{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook explores SNMR data and demonstrates an approach to estimating aquifer properties for a near surface aquifer. Presently it is a little messy but I intend on cleaning it up post AGC conference\n",
    "\n",
    "Neil Symington\n",
    "\n",
    "neil.symington@ga.gov.au"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sqlite3 import dbapi2 as sqlite\n",
    "import pandas as pd\n",
    "from shapely.geometry import Polygon, shape\n",
    "from shapely import wkt\n",
    "import numpy as np\n",
    "from hydrogeol_utils import spatial_functions, AEM_utils, plotting_utils, borehole_utils, SNMR_utils\n",
    "from geophys_utils._netcdf_point_utils import NetCDFPointUtils\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import math\n",
    "import sqlalchemy as db\n",
    "from sqlalchemy import create_engine, event\n",
    "import netCDF4\n",
    "import rasterio\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\",category =RuntimeWarning)\n",
    "from hydrogeol_utils.db_utils import makeCon, closeCon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function extracts the K profile using the SDR equation\n",
    "\n",
    "def SDR_K(df, C = 4000):\n",
    "    '''\n",
    "    df: dataframe containing GMR inversion data\n",
    "    N: empirical exponent for water content when estimating the water content\n",
    "    C: empirical constant for estimating water content\n",
    "    '''\n",
    "    return C * (df['Total_water_content'].values) * (df['T2*'].values)\n",
    "\n",
    "def objective_function(df, C, equation= 'SDR', loss_fn = 'L2'):\n",
    "    if equation == 'SDR':\n",
    "        # Find an estimate of log10 K\n",
    "        df['SNMR_K'] = SDR_K(df, C = C)\n",
    "    if equation == 'TC':\n",
    "        # Find an estimate of log10 K\n",
    "        df['SNMR_K'] = TC_K(df, C = C)\n",
    "    # Retrun the sum of the squares\n",
    "    if loss_fn == 'L2':\n",
    "        return ((df['SNMR_K'] - df['K'])**2).sum()\n",
    "    if loss_fn == 'L1':\n",
    "        return (np.abs(df['SNMR_K'] - df['K'])).sum()\n",
    "    \n",
    "def point_within_bounds(x,y, bounds):\n",
    "    if (bounds.left < x) & (bounds.right > x):\n",
    "        if (bounds.bottom < y) & (bounds.top > y):\n",
    "            return True\n",
    "    return False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Our next tak is to bring the SNMR data in \n",
    "\n",
    "# Connect to the SNMR database\n",
    "\n",
    "# Now we want to know the screened interval and other spatial information so we open up the \n",
    "# borehole database\n",
    "\n",
    "\n",
    "SPATIALITE_PATH = r'C:\\mod_spatialite-4.3.0a-win-amd64'\n",
    "\n",
    "# Add spatialite dll to path\n",
    "os.environ['PATH'] = SPATIALITE_PATH + ';' + os.environ['PATH']\n",
    "\n",
    "\n",
    "DB_PATH = r\"C:\\Users\\PCUser\\Desktop\\EK_data\\SNMR\\East_Kimberley_SNMR.sqlite\"\n",
    "\n",
    "engine = db.create_engine('sqlite:///' + DB_PATH, module=sqlite)\n",
    "\n",
    "@event.listens_for(engine, 'connect')\n",
    "def connect(dbapi_connection, connection_rec):\n",
    "    dbapi_connection.enable_load_extension(True)\n",
    "    dbapi_connection.execute('SELECT load_extension(\"mod_spatialite\")')\n",
    "\n",
    "\n",
    "connection = engine.connect()\n",
    "\n",
    "query = \"\"\"\n",
    "\n",
    "SELECT\n",
    "              s.site_id,\n",
    "              s.Field_ID,\n",
    "              s.mid_X,\n",
    "              s.mid_Y,\n",
    "              a.acquisition_id,\n",
    "              s.geometry,\n",
    "              s.elevation,\n",
    "              a.pulse_sequence,\n",
    "              a.pulse_length,\n",
    "              imm.Depth_of_Investigation\n",
    "              \n",
    "FROM \n",
    "\n",
    "             sites as s\n",
    "             JOIN acquisitions as a on s.site_id=a.site_id\n",
    "             JOIN inverse_model_metadata as imm on a.acquisition_id = imm.acquisition_id\n",
    "\n",
    "WHERE\n",
    "\n",
    "            a.pulse_sequence == \"FID\";\"\"\"\n",
    "\n",
    "df_acquisitions = pd.read_sql_query(query, connection, index_col = 'acquisition_id')\n",
    "\n",
    "# Since we are interested in the shallow resolution we will use the shorter pulse SNMR\n",
    "\n",
    "\n",
    "df_acquisitions = df_acquisitions.loc[SNMR_utils.choose_snmr_site_acquisition(df_acquisitions,\n",
    "                                 pulse_sequence_criteria = ['FID'],\n",
    "                                 pulse_length_criteria=\"max\")]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove some nulls\n",
    "\n",
    "df_acquisitions.at[:,'elevation'] = df_acquisitions['elevation'].replace('-9999.0', np.nan).astype(np.float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "40"
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "inv_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now extract the SNMR inversions\n",
    "\n",
    "acquisition_ids = df_acquisitions.index\n",
    "\n",
    "cols = ['Depth_from', 'Depth_to', 'Mobile_water_content', 'Bound_water_content',\n",
    "       'Total_water_content', 'T2*', 'acquisition_id', 'inversion_id']\n",
    "\n",
    "df_inversions= SNMR_utils.extract_snmr_inversions(acquisition_ids, connection,\n",
    "                                                  mask_below_doi = True)[cols]\n",
    "\n",
    "# And the inversion metadata\n",
    "\n",
    "# We also have an inversion metadata table, here I will pull out the metadata for \n",
    "# the inversion above\n",
    "\n",
    "# Find the site_id\n",
    "\n",
    "inv_id = df_inversions.inversion_id.unique()\n",
    "\n",
    "query = \"SELECT * FROM inverse_model_metadata;\"\n",
    "\n",
    "\n",
    "# A bug means I cannot pass params but I can pass (254,255)\n",
    "df_inversion_metadata = pd.read_sql_query(query, connection)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>site_id</th>\n",
       "      <th>Field_ID</th>\n",
       "      <th>mid_X</th>\n",
       "      <th>mid_Y</th>\n",
       "      <th>geometry</th>\n",
       "      <th>elevation</th>\n",
       "      <th>pulse_sequence</th>\n",
       "      <th>pulse_length</th>\n",
       "      <th>Depth_of_Investigation</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>acquisition_id</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>MillProf</td>\n",
       "      <td>466505.9100</td>\n",
       "      <td>8316696.170</td>\n",
       "      <td>POLYGON ((466296.2798914838 8316647.306548808,...</td>\n",
       "      <td>51.909279</td>\n",
       "      <td>FID</td>\n",
       "      <td>40.0</td>\n",
       "      <td>45.148358</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1</td>\n",
       "      <td>saltflats_st1</td>\n",
       "      <td>468448.6660</td>\n",
       "      <td>8354121.638</td>\n",
       "      <td>POLYGON ((468500.9411914838 8354169.254548808,...</td>\n",
       "      <td>3.142112</td>\n",
       "      <td>FID</td>\n",
       "      <td>40.0</td>\n",
       "      <td>24.773026</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2</td>\n",
       "      <td>2pm</td>\n",
       "      <td>468863.3010</td>\n",
       "      <td>8347878.272</td>\n",
       "      <td>POLYGON ((468914.1478914838 8347926.009548808,...</td>\n",
       "      <td>44.433735</td>\n",
       "      <td>FID</td>\n",
       "      <td>40.0</td>\n",
       "      <td>79.434384</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>3</td>\n",
       "      <td>7m</td>\n",
       "      <td>464216.0790</td>\n",
       "      <td>8350319.661</td>\n",
       "      <td>POLYGON ((464250.9278914838 8350216.810548807,...</td>\n",
       "      <td>18.368710</td>\n",
       "      <td>FID</td>\n",
       "      <td>40.0</td>\n",
       "      <td>79.434384</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>4</td>\n",
       "      <td>6m</td>\n",
       "      <td>458571.4590</td>\n",
       "      <td>8345568.043</td>\n",
       "      <td>POLYGON ((458622.3728914838 8345611.113548808,...</td>\n",
       "      <td>62.885246</td>\n",
       "      <td>FID</td>\n",
       "      <td>40.0</td>\n",
       "      <td>74.430066</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>415</th>\n",
       "      <td>261</td>\n",
       "      <td>OK18_70</td>\n",
       "      <td>499206.4563</td>\n",
       "      <td>8286306.341</td>\n",
       "      <td>POLYGON ((499258.881652 8286248.88157, 499148....</td>\n",
       "      <td>19.483648</td>\n",
       "      <td>FID</td>\n",
       "      <td>80.0</td>\n",
       "      <td>71.189093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>416</th>\n",
       "      <td>262</td>\n",
       "      <td>OK18_75</td>\n",
       "      <td>522306.2338</td>\n",
       "      <td>8301356.083</td>\n",
       "      <td>POLYGON ((522358.636074 8301298.60235, 522248....</td>\n",
       "      <td>11.542414</td>\n",
       "      <td>FID</td>\n",
       "      <td>80.0</td>\n",
       "      <td>66.368816</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>417</th>\n",
       "      <td>263</td>\n",
       "      <td>OK18_76</td>\n",
       "      <td>500099.6412</td>\n",
       "      <td>8295239.253</td>\n",
       "      <td>POLYGON ((500047.209849 8295296.70649, 500157....</td>\n",
       "      <td>19.070965</td>\n",
       "      <td>FID</td>\n",
       "      <td>80.0</td>\n",
       "      <td>71.189093</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>418</th>\n",
       "      <td>264</td>\n",
       "      <td>OK18_78</td>\n",
       "      <td>500110.0260</td>\n",
       "      <td>8293822.738</td>\n",
       "      <td>POLYGON ((500052.571133 8293770.30749, 500057....</td>\n",
       "      <td>19.967169</td>\n",
       "      <td>FID</td>\n",
       "      <td>80.0</td>\n",
       "      <td>61.771193</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>419</th>\n",
       "      <td>265</td>\n",
       "      <td>OK18_77</td>\n",
       "      <td>500111.1024</td>\n",
       "      <td>8295746.223</td>\n",
       "      <td>POLYGON ((500053.648477 8295693.79112, 500058....</td>\n",
       "      <td>18.596773</td>\n",
       "      <td>FID</td>\n",
       "      <td>80.0</td>\n",
       "      <td>66.368816</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>255 rows × 9 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                site_id       Field_ID        mid_X        mid_Y  \\\n",
       "acquisition_id                                                     \n",
       "0                     0       MillProf  466505.9100  8316696.170   \n",
       "2                     1  saltflats_st1  468448.6660  8354121.638   \n",
       "3                     2            2pm  468863.3010  8347878.272   \n",
       "5                     3             7m  464216.0790  8350319.661   \n",
       "7                     4             6m  458571.4590  8345568.043   \n",
       "...                 ...            ...          ...          ...   \n",
       "415                 261        OK18_70  499206.4563  8286306.341   \n",
       "416                 262        OK18_75  522306.2338  8301356.083   \n",
       "417                 263        OK18_76  500099.6412  8295239.253   \n",
       "418                 264        OK18_78  500110.0260  8293822.738   \n",
       "419                 265        OK18_77  500111.1024  8295746.223   \n",
       "\n",
       "                                                         geometry  elevation  \\\n",
       "acquisition_id                                                                 \n",
       "0               POLYGON ((466296.2798914838 8316647.306548808,...  51.909279   \n",
       "2               POLYGON ((468500.9411914838 8354169.254548808,...   3.142112   \n",
       "3               POLYGON ((468914.1478914838 8347926.009548808,...  44.433735   \n",
       "5               POLYGON ((464250.9278914838 8350216.810548807,...  18.368710   \n",
       "7               POLYGON ((458622.3728914838 8345611.113548808,...  62.885246   \n",
       "...                                                           ...        ...   \n",
       "415             POLYGON ((499258.881652 8286248.88157, 499148....  19.483648   \n",
       "416             POLYGON ((522358.636074 8301298.60235, 522248....  11.542414   \n",
       "417             POLYGON ((500047.209849 8295296.70649, 500157....  19.070965   \n",
       "418             POLYGON ((500052.571133 8293770.30749, 500057....  19.967169   \n",
       "419             POLYGON ((500053.648477 8295693.79112, 500058....  18.596773   \n",
       "\n",
       "               pulse_sequence  pulse_length  Depth_of_Investigation  \n",
       "acquisition_id                                                       \n",
       "0                         FID          40.0               45.148358  \n",
       "2                         FID          40.0               24.773026  \n",
       "3                         FID          40.0               79.434384  \n",
       "5                         FID          40.0               79.434384  \n",
       "7                         FID          40.0               74.430066  \n",
       "...                       ...           ...                     ...  \n",
       "415                       FID          80.0               71.189093  \n",
       "416                       FID          80.0               66.368816  \n",
       "417                       FID          80.0               71.189093  \n",
       "418                       FID          80.0               61.771193  \n",
       "419                       FID          80.0               66.368816  \n",
       "\n",
       "[255 rows x 9 columns]"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_acquisitions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAHgCAYAAAC7PehpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3Qc1dnH8e/d3nfVqyVZzU3uNm64YVpCM91UQ95AEhKIk9BLQiiBkEZICAmQBEgINSFAIKE7odlgm2LAxrjItlzVZclqu3vfP3YlbCzLK6Ft0vM5Z460s7Mzj433x8y9d+4orTVCCNFXhngXIIRIThIeQoh+kfAQQvSLhIcQol8kPIQQ/SLhIYToF1O8C4hEenq6LioqincZQgw5K1eurNFaZ/T0XlKER1FREStWrIh3GUIMOUqpzQd7Ty5bhBD9khRnHkIkkl27dlFXVxfvMqImNTWVrKysQ24n4SFEH9XV1VFeXo7RaIx3KQMuEAiwbt06CQ8hosVoNMJ/+9gON3dKdIoZQH0JRAkPIeLogQce4OOPP8ZgMDBx4kQWLVrUr/0sWbKEO++8E4ClS5fS0NDAwoULufHGG1myZAl/+9vfqK2txe/38+Mf/5j77ruP+vp6WltbmTdvHnPnzu3zMSU8hEgA27dv57jjjuOjjz7ir3/9KwDnnXce9913H3feeSd33nknCxcu5JJLLmHRokWsWbOGG264gSuuuIIxY8ZQX1+/3/6effZZKisrWbZsGQCXXHIJHR0dXHrppbS0tLBp0yZ+8pOfAHDZZZdJeAiRjM455xzGjh3L5ZdfTk5ODl//+tcBeOqpp7q3CQQCABQXF3P++eezZMkSVq9ezZQpU7jwwgu7Q6LLCSecwMKFC2loaACgra2N6667jmuvvRatNQMxFYeEhxBx9vDDD+NwOJgyZQrjx4/n/vvvB+D8889n2bJl/OY3v+GNN97g1FNPxWT6/Cs7duxYnnjiCR5++OFDhsGZZ57JmDFjePHFF7nooosoKirijjvuoLW1lVNPPbVfdatkmAxoypQpWgaJiUSxZs0aRo0aFe8yombfP59SaqXWuseWXhkkJoToFwkPIUS/SHgIIfolLuGhlDpWKfWpUmq9UurqeNQgxJelVN+WnjzwwAOcffbZANx1111885vfPGCbG2+8sbvXBOCxxx4DQmM7DmbfzyxZsoTm5mZuvfVWLr/8cp588kmCwSDXX389d955J1dffTXbt2/v858/5r0tSikjcDdwFFAFvKuUekZr/UmsaxEiEaSnp7Nx40bq6uqw2Wzs3LmTn/70p/h8PubMmQPAvffey7Zt2/jud7/L22+/zZlnnglATU0Nt956K+np6YwaNYpTTjmle7/33HMPdrudqqoqXC4X1113HZs2beIvf/kLLpeLmTNn8tWvfpU9e/Zw++23c+utt/ap7niceRwGrNdab9RadwCPAifFoQ4hEsKiRYu45JJLOOaYYwD43//+x3HHHcfVV1/dPdbjtNNO47zzzuPll1/e77NLly7FZDKRl5dHVVXVfu9961vfYsmSJeTn5wNQWVnJb3/7W6688soDunbVwU6NehGP8MgDtu7zuiq8bj9KqYuVUiuUUiuqq6tjVpwQseZwOPjVr37FjBkzAJgzZw7PPfcct99+OyeffDIATzzxBA899BALFizY77Nz586lvb2dlpYWJk6ceNBjNDY2cvLJJ5OTk8Nrr73G0UcfzRtvvMGvfvUrbr31Vi655JI+1x3zcR5KqdOBY7TWXw+/Pg84TGt96cE+I+M8RCKRcR4h8RhhWgUM2+d1PtBra01DbZB/PtTc60611nR2ttPW0UJHeytt7S20d7TS0bGXQKAFaCUYbCGoQ687O1sIBttI8bn45gUX0NnZSae/E78/0P17Z6c//DO83t9Jp98fWgJ+/IFA6PfOrnWh97+4LdDdYrb4ggsYPnx4///2hEgQ8QiPd4EypdRwYBuwCDi7tw9sqPyU0782DoNqBdUJ+FEqiNFgwGy2YLXYcdgdeFweXC4XVqsJm9WM3WrGZjdjtZixWC2YTCbMZjMWswWLzYrd5sJksfDg35/AbDZjMpsxW8yYzZZ9flowmU2YLRbMLjtOsxlzeOna3xeXfdebTCYMBukRH2wCgcCgnc8jUjEPD621Xyn1HeAFwAj8SWv9ce+fGoU/8A7QQOgkZQewC9iFxV9DZ2cNrW311Dc0olQbJmMHRqPGZAKjke6foXUKi8WEz+clJdWD1+siLc1DSqoXn9WC1eUiJSWFlJQUvF4vdrsdq9WKzWbDarVitVr3u79ADD2pqamsW7cu3mVETWpqakTbJcW9LUpN0TCQbR5thAJoJ6FAqseg6jGbGjCbGjAZ92Aw7MFoaMFg8GM0dGIw+DEYAijlx2gIhsPo81CyWszYbWbsNgsOhxmH3YLNZgmd5VgtWCxWLJbQmYxGYzAaQ6Fkt2O12bDarNhsod9tdhvW8Pqu0No3vPb9vT+t5EJEKtHaPBKADRgeXkKCGto7Q0vfaaAVqAc2Ap9hMq7Dbl2HUh+hdR0uu4PslDRK87L4+vFHMbGsnNrGRmqbGqnb00jt9hq272mitmtpCr23Z+9e2js7aetox2gy4nK6cDqduFwu3B4PqWlplI0o53vf//4A/L0IEbkhGh79sQV4BqetCqt5ByZTPRZzOyZTOyZjJ2ZTEKtV4fPYSfU5SfE6cLnScDqGYbKYQw2mhtBQw5V7qlm7rhWX240zxYWrIIMit5sKt7s7GFyuUEhYLJZ4/8GF6JGER8R24LCuwmjcTUDXogNtBNFYlQmDyYI2GggYFEGzHasvnazifIpLCykqKqK0tJRhw4ZJw6kYVCQ8IjaNve3ToL23bYKEeqI3ApswGbZgsSzDZNyJwdCAyaQxmzRmM5jNCrNZYzYrvD4PublZFA3PpaSskNLSUo444oiY/KmE6C8JjwFlAArCyzz8QfC3HeozQUINt+uBSoyGSiym3/Ls8786YDShEIlEwiPuDEBueJlDIAitHfm89MJrEh4ioclFeMJpw2H7DzXVjfEuRIheSXgkDD9262UMy5nB7+45kfv//Jt4FyREryQ84i6IzfJDcjInc9vtI9iy/T0WX3BuvIsS4pCkzSMqmoF1wHqU2oTNsgWTcTtGYx0mU3CfHhcwmzVnnnUcP7zxPenKFUlFwuOQurpf1wMbMZkqsZqrMBp2YjQ2YzaFQsBiobsL1uG0kJOdSeHwPMpGFFJSchwlJSUMHz5cBn2JQWMIhkcboSD4DNiIzbIVs2kbBkM1JlMAszmI2QQWiwqFgUWRkuIlLy+bkrJhlJWXU1x8LOXl5WRmZsb5zyJE/AyR8Lgeo/FxFLUEdQvoAMpgwGqxYbM5cTk9pKVm4nZ7cTrsOF0O3G4nHq8Tn8+NL9WD2+3uXoxGI9XV1fj9fnw+HzabTS45xJAzRO+q/aIgsJfQHbYNQGN42RNeGjEamjEa92AytGAwtKAMe1G0otRelOrEoMBg6Fp0+Kc6YJ3RqD5/z6gxGg047HbsDjsulxOXy4Hb48DrdeNNcePzhYLL4/F0L13TBfh8PgktEVVyV+0hGQBXeMnvcYtAMLR0DPixOwgF1b6h1RRe9qBUA2bjVgyG5vAUAS0oFQotaMeggvsE1IHhZTSG1imlw78rDOF1VqsVh92O0+nA6bLj9rjweF14vU5SUr24XC68Xm93ePl8PrxeL6mpqdhstgH/mxDJRcIj7ixARng5kNbQ4Y/GcYOEeoX2Da2u4AqdcZkMWzEaP8FoaMZg2Bs6yzK0hs+4vni21XNo7XcWFl5nMhlx2B3YHTaczs8vEb0+Nz6fC7fHTWZmZr8fwCxiQ8JjyDIAnvAC4CfUmNwaXvbgDzbiDzYATRgMjZiNTeHLtyYMqpmgYS86sJdgcC/QGQqU/QIDjAaFMoBBaQxGhdGoaKeTpsbdaO3HYrVgs9mx222h6QicDrw+J1m5aZx88slyWZbAJDwShubzL3DXl7htn6UZo6EBs6kRs7ERo2EPBkMTBtWCUl1tMB0o5UcZOlF0ggrN9br/Fznc9mIAg1F93hZjMmIxm7BYTJgtJux2Gx63E7fHidfrxO0JNRZ7fAV4vF48Hg9er7e7/cXpdGK32wflvJ6iZxIe+9FAJ/t/aff/EpuNezGb9mIy7sVkbOz+v7BSTSj2YjB0YjB0hr7EqjO8+AkFQyeoQKjRVH3+f2dD+LXJbMRsMmGxGLFYzVgsJqwWExarCYfDitPlwu1y4XK78Xg9uD15eH1ePOG2CIfDgd1u7/5pt9sxm83x+ssUg1yShIcf+ICuL3LoC9wa/hK3YjK2YjS2YTTsxWjowGTUGAydBIN14Wvzri9z1xc49FMT+lLve91uNhkxW0xYzCbMFmPoy2sxY7V2faHDs6NbLTgdrvCXOAOPpwS314PD6cThcnV/eb/4ZbZYLDLvqBgUkiQ82kh1/wu3oxKDYQtmUyseh4WcNC/jivOZVJ7PiGFFFOfmYbOEJgXeWVvDX157kQ83b2Lzrh24vV4mTp7EtJkzGTduHG63u3tmdLmuFqLvknycRzOwFbOxEp9rE057FV5nJx6nH68jgNcVpCQ3hfEl6RTlpNHW0UFlzS7W7dxOwGgAi5nUzAzGTZ7E2PHjZcSoEF/Q2ziPJA+PQwkSer7LVlz2jXidlTjtdXidge7F7fDjcQbAuJegyYDb58PssFM8opxxkycxatQoGdMghqwhHB6RaAW2YjRU4nNV4nJsxWNvxWKpQ+vd+HUdLreVrCwv+QXZzDj8cGbOnk1BQYG0XYhBT8LjS9FADVCJ1fwyBsM/6PSvIRhsIysrm2ee/SdTpvT4dytE0pPh6RHbDryFxbQSm/UTTKYabDawWcFmh+zsdCZNPoppM65i+vTp5Of3PJRdiKFgiIVHPfA2JsO72G0fYTLtwGbVoYCwQ2qqh/ETRjJj1mSmTfsaJSUl0hMjxEEMsvDYC7yDQb2Dw/4hJuMWbNZAKBxs4PHaGVNRzvQZE5g+42QqKiokHITopyQLjw5Cg8WW4bR9gMm0Aau1A7tNYbNpnC4L5eXDmTZzAtOmfZvJkyfLzF1CREmShMeHGAyhJ86PHDGGk085gcNnn81hhx2Gy+WKd3FCDElJ1NvyMPBfnPZlWEwbsNn82O0ah9NASWkhc+dPZf78eYwbN04uRYQYIIO8q9YPrMRg+C8O2wos5u2hULGDy21m7NgRHHHUTObNm0dhYWEsyxYi6Q3y8OjNXuB1zKbXsVvfx2Kpw24Hux1SU11MPWwc8xfMYu7cuaSmpg502UIkvSEcHr3ZDSzFZnkLq+UjbNbW7mDJzctgxqyJzD9iNjNnzpTh6WLIkvDos8+A13DalmE2b8Bu82O3g8OhuttXzlx0Brm5uTGsSYjYk/AYMH7gPQzqZSZPeJF3Vr0W74KEiCoZnj5gTMBUgnoEqJfjXYwQcSV9mv3SgUnm6hRDnIRHv3RgNMlJmxjaJDz6pVPOPMSQJ+HRLx0YjPJXJ4Y2OffulzxWrGglK30WdlvXoDM3k6dWMGfedObMmSPzoYpBT7pqB8x24A1sluVYzKuxWlu6g8XrszNh4mjmzJvG3LlzZRIhkTRknEfc1QKvYzEvw2ZZjcXSiN2msdvB5bYwdtxI5s4/jDlz5lBSUhLvYoXoJuGR0JqANzGb3sZu/QCzubY7WJwuE6NGlTFn3hRmz5nNqFGj5I5hEVMSHklrL7AMo+EtHPb3MZt2dgeLw2mkrKyI2XOnMHvO4UyYMEGCRQy4uISHUupPwPHAbq11RXhdKvAYUARUAmdoresPva+hGh696QDexWB4E6dtFSbTNuy2YOjmPoeiuHgYs+ZMZvbsWUydOhWTjEsR/RCv8JhD6JFuD+0THncAdVrr25VSVwMpWuurDr0vCY++8QPvY1Bv4LCtxGzejM0a6L5rOC8/izFjihlVMYLRo0dTVFREWlqa3D0sDhC3yxalVBHwr33C41NgntZ6h1IqB1iqtR5x6P1IeAycILAG+AiH9RNs1g1YTDWYzXuxmP1YrRq3y0J2ppeMDA8+nw+Pz0dGViZpGRmkZ2WRnpFBeno6aWlpMkfsIJdIN8Zlaa13AIQDRAZDxJwBGAOMYW877G3vaZtGYDM2y0Z8ro04bdtx2ddgs7yLxdyC1dJBms+Cy2XE7XLhdrsxmE1opcBoQBmNWGxW0jMzQ0tWJmnp6d2BI5dQg0PC/ldUSl0MXBx6VRDXWoYeLzCOto5x7Kzr6f02YAsm4yZ8ro24wg8Y97kCeF0dpLg0xblm7P5GTE3ttGzazs6Odmqam6jd00QAwGgAQ3gxGtBKYXc6usMmPSuL9HDgpKSkYJTbARJOrMNjl1IqZ5/Llt0H21BrfS9wL3RdtojEYQPK8QfKqWmEmsYvvh8AtqNUJV7nRjyOTbgdLXhdAXxOhdcVoCjbyfiSTEYWZFGYlYPP7WZvWxu1TQ3U1jdSs/lD3m9uomZPE3XNe7oDRxmN3cGjlcLpdpOeGb6cysrsDhyfzye9T1EW6/B4BlgM3B7++XSMjy9iwggMQ+thNDTPpqH5i+93Pf93My77RrzO5bjsdXhdAbzOAD5ngJw0K2OLM5laXE5hVjZZqWkHPFhca83etjZqGhuoqW6gZkMVlS17qNnTRH1LM9qgQpdS4bMbZTSiDQqPz0d6ZsbnbTjhwPF6vfLw8j6IZm/LI8A8IB3YBfwI+CfwOKHrkC3A6VrrHk+M99+XNJgOPXuAzVjNXZdGO/G6/PicAbyuAOkexaiidMaXZFCck0NeRibmCNpStNY0t+4NBU5jAzVNTdSEz3CaWveiDWqfSylD91mOLzXl80uqzEzS0tJIT0/H7XYP6sCRQWJiEGoHqjAYNpHi2ojbsRWPs707XFJcAUrzU5lYmsFho8rJTOn/7Phaa5paWj4PnD2N1OzZQ01zKHAwGFAmY3cbjjYolNFISlpqd+BkZGUxbty4pGsslvAQQ8weXPZnyEp5ltJ8I7d+/QQml5fGtIJgMEhD857uwNnV2MC7lRtwZKRx5uLzKCsvj2k9/SXhIYaAKjJ8/2B49meMLzVx0XGTmTJyRMJdUjS1NPPY66+xvr6GMVMnc+qiM3E6nfEu66AkPMQgpFHqPQoyn6E4p4ZZY71cdPwMCrKy411YxD7auJ4n33kTv93KiWecxtRp0xIu7CQ8xCDRgdX8CoXZr1CWt5eFhxdx+rzD8Cb5w847Ojt5ZvkbvLt5I/nlZZx1wfmkp6fHuyxAwkMktTp8rmcpzFrF6MIAFxw7lvmTxkfUs5KMtu7eySNvLKUu0MH8rx7LkcccE9cBchIeIslsICft7wzP3sqkcivfOGEqY4YXJ9wpfTQFg0GWfrCKlz75AE92FosuXMzw4cNjXoeEh0hwQYyGZRRkPU9JbgMLJqVzwbHTyU5LjFP3eKvf08Sjr79KZVM9E2ZO5+TTT4/ZHdASHiKBdI0urcJpW0dB1nLK89s4Y14ZJx0+BafdHu8CE9o7az/m7uefxpuZwVU/uoG8KM+Hm0h31YpBLUBoMHEVDttmvI7NOGzVeJwB3A4/HkcQrzNIYZaH0UXpjC7KYlzxmXLTW1hrexvbqqupqt5NVX0tVXU1tPo7wWQEoxFMRsw2K/OOO5b8wkJcbndc65XwEBHqJDRD/Fbcjs14HFtx2OrwOIN4HAE8jgA+l6Y4N4WKojRK8jLIz5hAqkfuFwFoaW2lqnoXVTXVVNXVsq2+hvZgoDsUtMGA3eUkv7CQ/HFlTBt2BKfk5eFwOOJd+kFJeAhCt9hvQ6mteJ2bcdm34LLvwe0I4HEE8TgDpLoVZfmpjClKY3hOBvkZ0/E4nRIMhAZ+VVXvDgVDfQ3b6+vo1Dp0X0z4rMHpcYeCYdJoZg8bRm5ubtLP3CbhMeg1A1UYDVvwujbjslXhcrSGQsHhx+MMkOYxMqIgnYqiNAqzM8lLnydtD4TuaWlo3tMdDFvratjRUEdQqX0uJQx4UlLILywgf9p4FgwbRk5OzpCYYU3CI2lpQjN+VWE2bcHn3IzDVoXb0YHHGT5jcATISrUysiCdMUVpFGRlkps2BusQ+Id9KFprahob9jljqGVXYz3aaOi+lMBoJCU9LXTGMHsqY/Pzyc7OTrqb26JF/hYS0uc9ElbzZryuzTitO/E4/eFLiQAeZ4D8DBcjC9IYXZRGfsYwslMnDtrBU30RDAbZXV8XCobaGqrqa6je09R9ptAVDumZWQwbXkj+hMOZmp9PZmamTCDUB/IvLea6eiS29dAjEQoGrzNIQZabMYXpjChIJz9jBJkpM+QfNhAIBNhZVxsKhroaqupqqGtp3q9HwmA2kZmdTX5RIYWHjWVWfj7p6enSPjPAJDwGVCewA6jCba/E7diCs6tHwtkVDOEeieFplOZlkJ8xXnokwjr9fnbU1oR6JepCXZWNbXtDoRAOBqPFQnZeLvmFBZTNmsT8/HxSUlLk7y8OJDwGgMd5HeX57zGupIzSvFTGDpceiZ50tTN8VrWF9bt2UFm9m4BBgTkUDiabldz8fPKLCxkzbzrH5Ofj9XrjXbY4CAmPAdDU8mOqqu9iUtlWlpy2ALs1ubvgvgytNbvqalm/rYr1u3ewuaYrIExgMpKZk0PpyJHMPGYuZxUWYjab412y6CcZnj6gKhk7/DZu+tpEFs7ucUTvoKC1ZkdtDeu3bWX97h1sqa0hGB7ToMxmsnJzKB01ktLycgoKCqR3IonJvS0xpfG5HmfOuFe5/8rTyPClxLugfgkGg2yvqe4+g6iqryVoNKDMJjCZyBmWT+nIEZSNGEF+fr4MMR+kJDziopHSvFv45olevn/GMQnZ7hEMBqmq3s367VV8tnM72xvqui8vlNlMbsEwysJnEHl5edLbMwRJeMSRzfI6U0b8iXt/cDyjigpjfvxAIMDW3bv4bNtW1u/eyc7GejCb0EYjBquZ/MJCykaPorSsjJycHAkIsR8Jj7jrIC/jl5w6u4Y7vnnKgI/w9Pv9bN61k/XbQ5cYu5oau88gjFYrBcVFlI4cSWlZGdnZ2Ql5FiQSk4RHgjCozxhXcgc//cZ0jp46vk+f7fT7qdy5nfXbt7F+947QiMlwQJhsNgpLhncHRGZmpgSEGBASHglFk+p5kAWTlvOHH5xOitvT/U5HZyebdmxj/Y5tfLZzB7V796DM5tA8Dg47RSUloUbK8nLS0g58/KIQA03CIyHVMCzzXE49Kgt3SgqYTFicdoaXllI6aiRlZWWkpvb/KWdCDASZSSzhdJKZ8n2KSwJc/+tfkpaWFu+ChOgzaVqPMcVH5GZM44orK3jtzRclOETSkjOPGPI6byYv92mefu4xSsvK4l2OEF+KnHnERD0ZKXM4/bRtrF77jgSHGBTkzCPKrKYnyUi/kb88fBfzjjgi3uUIMWAkPKLGT4r7HCaM38NzL76LXeYEFYOMXLZExQdkpEzittsX8Orrz0twiEFJzjwGmMN2HQXDXuaVpf8hNzc33uUIETVy5jFgakhxz+SCxW2sWbdcgkMMehIeA8RsGk/esA7GjC2luro63uUIEXUyPH1AfYzJ+ARO++s47HtxuaGkJI/TzvgKp552Kj6fL94FCtEncm9LXK3CYvo7dttbOB1tuNwwYkQRZyw6joUnL8TlcsW7QCEOSsIjoQSBd7GYn8RhexenswOXS1NRUcaZZ5/A8ccfL70zImFIeCS8IPAWNsvfsdtW4nT6cbthwoRRLDrnJI499tgh8exTkXgkPJKSH/gvdus/sFk/xOUK4HErpkwdy6JzFnLkkUfKrOQi6iQ8Bg0/8CIO6z+x2T7G5dR4PAZmzJzIonNOZt68eTIHqRhQEh6DWhvwb5y2Z7BZ1+Fyabw+E7PnTOWsc05hxgx5xq3oPwmPIWcv8Cwu+7NYrRtxuyAlxcK8I6Zz9rmnMWnSJAkUEREJDwE0AU/hcvwbm3ULbhekptk46pjDOefcM6ioqIh3gSIBSXiIg6gD/o7b8W/stp24XJCR4eTY4+ZyzrlnUibzjgx5cQkPpdQw4CEgm1Bf5L1a618rpVKBx4AioBI4Q2td3/u+JDxiZzeKJ3G7bqVxz7Z4FyPirLfwiOaFrx/4gdZ6FDAd+LZSajRwNfCK1roMeCX8WiSMTDTzGTtmXLwLEQkuauGhtd6htV4V/n0PsAbIA04CHgxv9iCwMFo1iP5xO67jF7++Md5liAQXkyZ3pVQRMBFYDmRprXdAKGCAzFjUICLn825h2rRp8S5DJLioh4dSygX8HViitW7qw+cuVkqtUEqtALnFPVYUf2fuvInxLkMkgaj2tiilzMC/gBe01r8Mr/sUmKe13qGUygGWaq1H9L4faTCNlRTPHNZ+9iSZmXJCKOLUYKpCD1L9I7CmKzjCngEWh39fDDwdrRpEX7WRntYuwSEiEs07q2YB5wGrlVLvh9ddC9wOPK6U+j9gC3B6FGsQfWAx3cbF3zwr3mWIJCGDxES3jNSpbN/1ttytK7rFa5yHSCqbGF7kleAQEZPwEAC47Ndw28+ujXcZIonI/2aGrDrgVWyW17FZVpOS0sAR8jhM0QcSHoNaEPgUxcs47cuwmCuxOzQOO/hS7Ew9bBxHHj2P+fN/LDO7iz6T8BgU2oA3MZtew25bhdXSgMMOdgfk52cxd/5hHHX0ZUyePFnaNMSAkX9JSWUn8DIO6+tYLGux2zqx28HlNlFRUc6Rx8xiwYKLKSgoiHehYgiQ8Eg4QeB9DOoVnI7lWEzbsTvA4YD0NA8zZk1kwVELmTt3Lg6HI97FiiFMwiNumoGlWMxLsVvfx2bdi90ODqdi+PBhzD9yGscccyOjR4+WKQNFQpLwiLpNwEs4bW9hsazHYQ9it4PHa2XixNEcdewc5s+/UoaEi6Qj4TEg/MAyTMZXcdhWYLHU4LCDw2kgOyuFw+dO4aijv8b06dPl4U1i0JDw6JM64GXs1texmD/Cbu/oDokRI4az4KiZHHnUL2TuTzEkSHgcIAisQfEKLsfbmM1bsduDOOyQkupg6mHjOeqYI5k//1Y8Hk+8ixUiboZweLQBr2MxLcVuXYXV2hhqsHQohhWExrl7R6UAACAASURBVEYcfcwP5BknQhzEEAiP7cDLOGyvYzF/it3mx+EAt9tMxbgRHHXM4RxxxLfIz8+Pd6FCJJVBFx4O20+x254Ot0Uo0tM9zJo1iSOPOYPZs2djs9niXaIQg8KgCg+b+Q5mzljGCy+/IZcaQkTZoAkPk/FhRo96npdeXRrvUoQYEgZFeCheobjwVyxfuSzepQgxZAyC8PiQvJwlvLd6udwxKkQMJfm3bQuZ6YtY+f5SuUlMiBhL4lbFBtJ8X+GNt56W+0KEiIMkDY8OUjzz+Nfzf5Sh4ELESRKGRxCvay5/euAmps+YHu9ihBiyki48vK7juPnWC1h48onxLkWIIS2pwsPtOJeLvzGVSy/7RrxLEWLIS5rwcNq+z0knObnj5zfFuxQhBEnSVavUTmbOrOQvf/tHvEsRQoQlxZmHYifnX3Ay1dXV8S5FCBGWJA+6HqGd1hOx21fhsO/F41YUFqQxYeJoZs2dw+iKCoYNGyY3wwkxwHp70HWShMcUDSv2WeMH1uJzv4TH/i52ezVWazPFBekUFuWQkZ1N+ejRjBpbQVlZGVarNV6lC5HUBmF49KQDWE2G723SPZ/hdVVjtzWSlWkmJzsdl9eDwWKhsLSYkRUVjBo9Wh6xKMQhDJHw6EkrSr1PVsrbZPq2kZ3aTqavjdJhNow2A50AZhOYTaRlZTJybAWjxowhPz8fpdQA/ymESD5DODx6sgejcQVZvmVkplSTk9bByAI7h4/NwWWzsKW+hqr6WoImI5hM2D2u7kug0tJSeXSCGFIkPA6pDrPpXbJSlpOVUk9OWidji90cN72MkQWFVFXvYs22LazftYNOBZhNGKwWikpLGFlRwchRo/B6vVGsT4j4kPDol11YzcvJTl1BVmoTuWmdTCpL4yvTyhhXUorRYGDzrh2s2bqZtTu20dTeBmYT2mQkMzen+xIoNzdXLoFE0pLwGBAaqMJhW05Wyntkp7aQm9bJtNFZfOWwckYVDsdoNKK1pqaxgTVbKlmzbSvbG+vRJiOYTTg8bkZUjGFURQUlJSWYzeY4/5mE6J2ER9RoYCMe5zIyfavJTm0jL93P7HG5HDVlBKV5+489aWltZV3VFtZs28KG3Tvx73MJVFxezsixFYwcORK32x23P5EQ+5LwiKkAsI4U19tk+NaQm95JfkaQBZOGMW/CCAqzcw64jPH7/VTu3MGaqs2s3VFFc0cH2hxqsM3Oz+u+BMrOzpZLIBFTEh5x1wl8Qrr3LTJ8G8hJ66AwS3H0lCLmjB9JbnpGj5/SWrOrrpa1WzezZnsVO5saui+BXD4vI8eOZeSY0RQXF8v8rSIqJDwSUhvwIZm+t8hMqSInrZ3hOUa+Oq2EWRWjSD/EALbmvXv5dOtm1mzbysbqnQQMCsxmjDYrJSNCl0AjRozA5XLF5o8jBiUJj6TRgkGtIiv1bTJTdpGb1klZnoWvTi9l+uiReCMIgk6/n43bq1hbtZW1O6po8XeCOXS2kltQELoEGj2azMxMuQQShyThkdQaMRnfJTt1OenenZw618f1553U571ordlRWxO+BNrKrj2NYDZRUjGGxRdfFIW6xWAg4TFIFGb9kBd/PonyYQP3UO7f//sZppx8PFOmHTZg+xSDR2/hIfewJwmn7QW+f4ZjQIMD4OJjjuehu++htbV1QPcrBr+IwkMpNU4pdaJS6pSuJdqFiX1Vc+Tkf3DZqUcO+J4NBgNXnHAqP7vplgHftxjcDtm/p5T6EzAO+BgIhldrQOYEjAlNxfDreejas6J2hGGZ2RRbnbz28svMP3LgA0oMTpEMDpiutR7d1x0rpWzA/wBr+DhPaq1/pJRKBR4DioBK4AytdX1f9z9UZKfexe+WHI7HGd0u13PmHcXlD/6ByYcdhsfjieqxxOAQyWXL20qpPocH0A4cobUeD0wAjlVKTQeuBl7RWpcBr4Rfix6YTas4/+gdzB4/KurHUkpx9cIzue2GH0X9WGJwiCQ8HiQUIJ8qpT5USq1WSn14qA/pkObwS3N40cBJ4X127XthP+oeAlqYMfo33HZx7JqXMnwpzMgt5Okn/x6zY4rkFUl4/Ak4DzgWOAE4PvzzkJRSRqXU+8Bu4CWt9XIgS2u9AyD8U55S3YPSvB/y6I/OjPmkzidOP5xl/3lJZqoXhxTJv8wtWutntNabtNabu5ZIdq61DmitJwD5wGFKqYpIC1NKXayUWqGUWgFD6x+yz/UoP76wiJy09Lgc/5pTFnHb9T8kGcYAifiJJDzWKqX+ppQ6q79dtVrrBmApobOXXUqpHIDwz90H+cy9WuspoQEqPd84NhgpVckJM9/g7CNnxK0Gj9PFiRWTePiBB+JWg0h8kYSHnVDj59GELle6Ll16pZTKUEr5wr/bgSOBtcAzwOLwZouBp/te9mDlZ2LpTdx3+dnxLoR54yay4Z1VbN2yJd6liAQVteHpSqlxhBpEjYRC6nGt9U1KqTTgcaAA2AKcrrWu631fQ2N4+rDMm3nu9tGMLS6MdykAtLa3cdWjD3Dn/ffKA7WGKLm3JQnYra9x4+LXufLsr8a7lP2s/Gwt7+6t45uXXRbvUkQcyL0tCa+O+RMe5oqzvhLvQg4wuWwkezZuYe2aNfEuRSQYCY8EYLN8jNelaOtoj3cpPVpywmk8fN8f412GSDCR3NsyBZgN5AKtwEfAy4dqpxCRa+uYzSOvDOPjypu569I5zJ0wJt4l7cdsMmGMdxEi4Rz0zEMpdYFSahVwDaEel08JdaseDryklHpQKVUQmzKHgiI+3HA/Z9zYwjd+/jDtHR3xLmh/wcRvGxOx1duZhxOYpbXucaIHpdQEoIxQj4kYEIrdDZdw33MbeWftrdy9ZD4zK0bGu6iQYPDQ24gh5aBnHlrru4EOpdT3DvL++1rrV6JW2RCmdTHvr7+PU26o49t3PkJHZ2e8S5LwEAfotcFUax0gdCObiDkDu+ov4/dPf5MZl9zLO2vWxbccuWwRXxBJb8ubSqnfKqVmK6UmdS1Rr0wAENRlrPrsfk66biffvesxOv3+uNShtCYoZx9iH5FMBjQz/POmfdZp4IiBL0f0zMjOuu9z9z/X8OZHd3Df5ccysbwkphW4bXaam5tloiDR7ZDhobWeH4tCxKEFgqNYue5+jr/mlyw6YiU//cYpMXtSnM/ppKGhQcJDdDvkZYtSKksp9Uel1L/Dr0crpf4v+qWJnhnZXnsFd/3jfGZ+5/es3rApJkf1OZzU18tskeJzkbR5PAC8QGiQGMA6YEm0ChKR8QcqeHft/Xzlqg1c9YcnCQQCUT2ez+6kQcJD7COS8EjXWj9OeOZ0rbWf0KPgRdyZ2FZzNb984ixmfece1lRGNEdTv6S4PTTUyaBi8blIwqMlfBu9BghPYtwY1apEn/j941m+5j6OvmIt19//VFR6RXwuF/W1tQO+X5G8IgmP7xOawKdEKfUm8BBwaVSrEv1gpqr6Ou549BQOv/Qe1m2tGtC9+1xuGurkskV87pDhobVeBcwl1GX7DWCM1vqQs6eL+Oj0T+btj//AkT/4gJsefGbAzkI8DieNDQ0Dsi8xOETS22IDLgNuBn4MfDu8TiQsC1t3/4if/PU45n73Hmobv/yX3mg0EkiEYfIiYURy2fIQMAb4DfBbYDTwl2gWJQZGe+c03lh9Bzf86d8Ds0O/tJOLz0UywmhE+KlvXV5TSn0QrYLEQMtg6fua9o4OrBbLl9uVhIfYRyRnHu+Fe1gAUEpNA96MXklioK3bupi7//nlb4BW/oDc3yK6RRIe04C3lFKVSqlK4G1gbqSPnRTxFwiO55FXtn/phzjl+lLZsWPHAFUlkl0kly3HRr0KEXVrthzLf955j69M6/8N0WXZuXy2bh15eXkDWJlIVpF01W4G0gjN63EikNbXx06K+GtpPYFfPL7yS+2jNDef9WvWDlBFItlF0lX7Q0IPb0oD0oE/K6Wuj3ZhYqAZWL1hFOur+j94LD8jk6rN8v8LERJJm8dZwFSt9Y+01j8CpgPnRLcsEQ27G77Oj/7c/4ZTg8FAsF3GeoiQSMKjEth3UJgV2BCVakSUuVi+xkNTS3P/dxGnmcxE4okkPNqBj5VSDyil/kzouS3NSqm7lFJ3Rbc8MdA2bL+Y2//2Qr8/r/yBL91rIwaHSHpbngovXZZGpxQRG4U893YzN38tgNHY90c55fhS2LFjB7m5uYfeWAxqkUxD+GAsChGxs3bLGfz1pTdZfOycPn+2LCuX9Z99JuEhIupt2aSU2vjFJRbFiejo8B/Off/q36McpLtWdInksmXKPr/bgNOB1OiUI2JD8XHlDN5Z8ymHjRrRp08Oy8xi63tvRakukUwiGSRWu8+yTWt9J/LYhaTX0HwWt/yl77coSXet6HLIM48vPODJQOhMxB21ikSMWFi1Lp+dtTVkp6X37aPSXSuI7LLlF/v87gc2AWdEpxwRS9tqTuGx157ku6f17falgISHQB76NKQZDNsYU9S3s45AIIByyERyopc2D6XUuUqp3t4vUUodHp2yRCykudcxoiC/T595+5OPmDFvbpQqEsmktzOPNEITAa0EVgLVhHpbSglNiFwDXB31CkXUuJ2byUuf3KfPvPrJB1z59TuiVJFIJgcND631r5VSvyXUszILGAe0AmuA87TWW2JTooiWVLcfgyGSOxQ+12ZU2Gxy2SIO0eahtQ4AL4UXMcikuPvW8LlxexXFo0dGqRqRbPr2vx0xiGhS+xgez61cznELF0apHpFsJDyGrGpKcj19+sTOvc3k5OREqR6RbCQ8hqz1TBsVeRA0NjfjzuzjYDIxqEUywtQKnAoU7bu91vqm6JUloi3F9SmjiiKfyPiFVe9wzIknRLEikWwiOfN4mtDkx36gZZ9FJDGvayOFWZGfeXywbTMTJk6MYkUi2UQyPD1fay2PXxhkfK5WLGZzRNsGAgEMTjtKqShXJZJJJGcebymlxka9EhFTqZ7IHx0po0pFTw565qGUWg3o8DYXhicAagcUoLXW42JTooiGvnTTyqhS0ZPeLluOH4gDKKWMwApgm9b6eKVUKvAYoQbYSuAMrXX9QBxLRKqR/Ax7xFvLqFLRk4NetuzzRLhb9n1CXNe6Phzju4SGtHe5GnhFa10GvILcHxMHG5g6MiuyLbfJqFLRs0jaPMbs+yJ8JhHR3VRKqXzgOOD+fVafROgJdIR/ypDFGHM7PmVccWTdtM+tklGlome93XJ/jVJqDzBOKdWklNoTfr2bUPdtJO4ErgSC+6zL0lrvAAj/zDzI8S9WSq1QSq0I3dArBkqKay3FuYe+FV9rzeaGOhlVKnrU22XLbVprN/AzrbVHa+0OL2la62sOtWOl1PHAbq11v56urLW+V2s9RWs9BTL6swtxEMNzqnFE0IZx93NPcfqF58egIpGMIhnnca1S6hTgcEK9L69rrf8ZwedmAScqpb5KaB4Qj1Lqr8AupVSO1nqHUiqH0JmMiBGD+oBT5x76rOP1jz7AUpjH9FmzYlCVSEaRtHncDXwTWE3oUZPfVErdfagPaa2v0Vrna62LgEXAq1rrc4FngMXhzRYT+SWQGADFuY/yta/M7nWbHbU1PPPJ+1z07UtiVJVIRpGcecwFKnT4AaVKqQcJBUl/3Q48rpT6P2ALoefAiJjopGJ4E077wbtpO/1+bn3qUX76+7tlRKnoVSTh8SlQAGwOvx4GfNiXg2itlxJ+xq3WuhZY0JfPi4HhtD3Ptxf2fn/K7X//G5ddfw1OpzNGVYlkFcllSxqwRim1VCm1FPgEyFBKPaOUeiaq1YkBVZr3Mgsmjz/o+/94639UzJ9N+Yi+PUVODE2RnHn8MOpViBioZVaF5aCXIp9u3czHzXXccPr3Y1yXSFaRPLflv0qpQqBMa/2yUsoOmLTWe6Jfnhgomb6HuPKseT2+19Laym9efo477783tkWJpHbIyxal1EXAk8AfwqvygUi6akUCGVnwGYXZBw720lpz42MPcf3tP8FkiuREVIiQSNo8vk1ozEYTgNb6Mw4yKlQkJqVWc8qc3B7fu/c/z3Li4nPIzs6OcVUi2UUSHu1a646uF0opE6HBYiJJFOc8wtePO3A+jmVrPsafmcLsefNiX5RIepGEx3+VUtcCdqXUUcATwLPRLUsMHD9jhjccMLajuqGeJz5YziXfWxKnukSyiyQ8riZ0Z9pq4BvA88D10SxKDByH9Xm+c/L+YzsCgQA3PfkwN97xUxkIJvotkt6WoFLqn8A/tdZye2uSKct/mQWTzt1v3c+eeoxvXXU5brc7TlWJwaC3W/KVUupGpVQNsBb4VClVrZSScR9Jo44ZY0z7PY/2X++8xfBpkxldURHHusRg0NtlyxJCvSxTw7fhpwLTgFlKqe/FpDrxpWSm/JUrz/q8oXTj9irerd7OmeeeE8eqxGDRW3icD5yltd7UtUJrvRE4N/yeSHAj8tcyPCc0Y1hrexu/+Pc/ue4WeVaXGBi9hYdZa13zxZXhdo/IHvgh4kapj1g4+/NBYTc9/leuufVmLBZLHKsSg0lv4dHRz/dEAjAbm3DaPg+K+tYWsmU6QTGAeguP8eG5S7+47AHkIVAJrsM/kz88u4lAIPRwpyXHLuQXt/4kzlWJwaS3OUyN4blLv7i4tdZy2ZIEPt60mN8/8yoAIwsKyfIrlr7ySpyrEoNFJIPERJLq8M/gvn9twe8PPR1u8YJj+ddfH6G2tjbOlYnBQMJjkPtk89f4zVMvA6CU4obTzuXmq68lPKukEP0m4THIdfqn8sC/t9MZPvvwulwsmjyTe397yDmsheiVhMcQ8Mnmi/nV4y90v54+agx7N23lg/fei2NVItlJeAwB/sAE/vpyNR2dnd3rvnvCqdz/q1/T0tISx8pEMpPwGCI+qfwGdzzyn+7XBoOBa046k1uvvyGOVYlkJuExRASCY3n01Xra2tu71+WmZ3B43nAef/hvcaxMJCsJjyFkzeZL+MnD/95v3VenzmDNG2+zaePGOFUlkpWExxAS1CN58r/N7G1r22/91aecxc9vvInOfdpEhDgUCY8hZu2W73Dzg8/tt85qsfDdY07kZzffEqeqRDKS8BhitC7lqTfaad67d7/15fkFDFMWXnnxxThVJpKNhMcQ9OnWy7jxgecOWH/u/KP5zyNPUF0ts02KQ5PwGJKKePatAE0tzfutVUpxw+nncMs118nwdXFIEh5D1LqqJdzwx38dsN7jdHHO1MO5585fx6EqkUwkPIasfJ5fbqRhz4GPHD5s5Gg6q3ayasWKONQlkoWExxC2ftsSrrm/5+d3XXr8KTxw1900Nzf3+L4QEh5DWg7/XpbJX15844B3DAYD1y48kysu+Q6Vmzb18Fkx1KlkaBhTaooGOYWOFrf9aWZWPMufrjqV3PSM/d5rbW/jt8//k4DPzaVXXI7T6YxTlSIelFIrtdZTenxPwkOENDM85zYu/IqJ6849br8HRQFs2bWT377wLJPnz+GMc86Wx1QOERIeImIW07tMKv8dv/vesUwsKzng/f+tfp+/r1rGhd+5hAmTJsWhQhFLEh6ij/xkp97FibO2cOd3TsNute33biAQ4E8vP8/mthaWXHs16enpcapTRJuEh+inTYwtvo1b/m8yJ86afMC7dU2N/Pq5p8geVc5F374Ek+mQz00XSUbCQ3wJmhT3I8wd/1/uu/x00n2+A7b4aNMG7lv6AieevYgFRx8dhxpFtEh4iAFQT2neLXzn5DQuO/WoAxpMtdb8483/8tbWjXznyssZXlwcpzrFQJLwEAPGblnK1JEPcd8VJ1A+bNgB70vX7uAi4SEGWDv5GT/nzPlN3HbxyZh7aOvYsmsnv3nhGSbPm8OZ554jXbtJSsJDRIXRsIbxJb/gF5fMZt7EMT1u87/V7/Pkyre58DuXMHHygY2uIrFJeIgoCpLuvZ+jp3zAPd8/HY/TdcAW0rWbvCQ8RAzsYmTBLVx1djEXHDu7xy2kazf5SHiImHE5nmPGqKf401Unk5+Z1eM20rWbPOIWHkqpSmAPEAD8WuspSqlU4DGgCKgEztBa1/e+HwmP5NJCUfbtLD5Gc8P5J2A0Gg/YQmvNU2/9jze3bJCu3QQW7/CYorWu2WfdHUCd1vp2pdTVQIrW+qre9yPhkYzMplVMLPsNv1tyDJNHlPa4jXTtJrZEC49PgXla6x1KqRxgqdZ6RO/7kfBIXgGyUn/DiTMqufPS03DYbD1uJV27iSme4bEJqAc08Aet9b1KqQattW+fbeq11ik9fPZi4OLQq4LJsDlqdYpY2EzF8Nu46cLxnDxn6kG3kq7dxBLP8MjVWm9XSmUCLwGXAs9EEh7770fOPAYHjc/1GLPHvcYfrzyNDF/P/9mlazdxJERvi1LqRqAZuAi5bBniGinNu4VvneTle6cfc9BLlK6u3ayRZVz07Uswm80xrlP0Fh5Rm8NUKeVUSrm7fgeOBj4CngEWhzdbDDwdrRpEovKyftvPuO7+ucy57HdU7e75IVOpHi8/PusC5viyufyib/DW/16PcZ2iN9GcADkLeEMp9QHwDvCc1vo/wO3AUUqpz4Cjwq/FENTWMZs3Vv+SOx59pdftKoaXcOcF3+KlRx5nzccfx6g6cShRCw+t9Uat9fjwMkZrfWt4fa3WeoHWuiz8sy5aNYhkkMLHlR2H3EopxfWnn8s9P/05dXXyTyYRyKMXRNxt3WXH7/cfcjuj0chNZ5zPD79/RUTbi+iS8BBxt61mOu9++mlE2/rcbi476jhuvvb6KFclDkXCQ8Td3vZ5PPrKRxFvX55fwKycAv78h3ujWJU4FAkPkQBS+LiyvU+fOHrSYbRu3ML/XnstSjWJQ5HwEAlhy+7I2j329a2vnMS/H35UHocZJxIeIiFsq57GynWRtXt0UUpx45mL+fmPfiwP5I4DCQ+REPa2z+PRV/s+hsNqsfDDU87mhh9cTjAYjEJl4mAkPESCSGX1xrZ+fTIzJZXF0+bws5tvGeCaRG8kPETC2NqPdo8uE0rKKbO4eOJvjwxwVeJgJDxEwqiqOYxVn63r9+dPmTmHyndWsmqF3EQZCxIeImHsbZvPo69GPt6jJz846Qz++tt72Llz5wBVJQ5GwkMkkP63e3QxGAzcvOgCbr7yatrb+zZ2RPSNhIdIKFt2WQkEAl9qH067nauOP40fXXEVyfB0gGQl4SESyraaaV+q3aNLQVY2J4wcx92//NUAVCV6IuEhEkrLALR7dJk1Zhzuxr3851/PDcj+xP4kPESCSWP1hr0DtrfFC45l+XP/Ye2aNQO2TxEi4SESzuZdti/d7rGv6047h7t/8lPq63t9tpjoIwkPkXC21RzGe599NmD7M5lM3HTm+fzwB5fLJEIDSMJDJJyWtvk89trAtHt0SXF7+M4RX+WW628Y0P0OZRIeIgGls35by4DvdcSwQopMdj744IMB3/dQJOEhElAQtyM6j5ucOLyU9Wv7duu/6JmEh0hA1RRle6Ky58KsHDZv3BiVfQ81Eh4iAVUxujAtKnv2ud00yqMbBoSEh0g4NssWyodlRO8A/oHrBh7KJDxEwvG6NpOfkRm9A0h4DAgJD5FwnLZdpHt90TuAhMeAkPAQCcfj8GMwRO+fptdmp6GhIWr7HyokPETCcTuie2ZQmJ7J5s2bo3qMoUDCQyQcT5TDoygji80b5VkvX5aEh0gwGo8zuo9QKMzKoXL9hqgeYyiQ8BAJpoG8dGdUj5Dm9VK7e3dUjzEUSHiIBBO9AWJdlFLS4zIAJDxEQjGbtjCqMD36B5Jb8780CQ+RUHzOLdEdINZlACcbGqokPERCcdq3kZ0a3csWAKfZSkvLwN/2P5RIeIiE4ra3YzKZon6cwrQMGevxJUl4iITidsbmcqIwPVPGenxJEh4ioUR7gFiXwuwcKjesj8mxBisJD5FQPM7YPOEtKyWVXdt3xORYg5WEh0ggzWT6zDE5ksFgQHdKd+2XIeEhEsg2RsZijEcXGSj2pUh4iIRhMGxlVEFq7A4o4fGlSHiIhBEM+thdH7uxFyk2B7W1tTE73mAj4SESyAje/GhbzI42aXgJ761cGbPjDTYSHiKBuKna3RGzo00oKeO95e/G7HiDjYSHSCg1TdEfXdrF43SxRx7D0G9RDQ+llE8p9aRSaq1Sao1SaoZSKlUp9ZJS6rPwz5Ro1iCSS21j7MIDgI7O2B5vEIn2mcevgf9orUcC44E1wNXAK1rrMuCV8GshAGhozqImhpMTu80WmpqaYna8wSRq4aGU8gBzgD8CaK07tNYNwEnAg+HNHgQWRqsGkXxqmir4ZHNlzI43saCY91atitnxBpNonnkUA9XAn5VS7yml7ldKOYEsrfUOgPDPGEzeIJJFp38Ub31UFbPjTSwtZ9Wy5TE73mASzfAwAZOAe7TWE4EW+nCJopS6WCm1Qim1IpRBYmgYxupNsRt7keb1UbdL5jPtj2iGRxVQpbXuivUnCYXJLqVUDkD4Z4//5bTW92qtp2itp0AUn1sqEoyB6gZjbA8pjab9ErXw0FrvBLYqpUaEVy0APgGeARaH1y0Gno5WDSI51TbGNjycRpPMKtYP0e4XuxR4WCllATYCFxIKrMeVUv8HbAFOj3INIsnU7bHT0dmJxRybO2zHDxvOhx98wIyZM2NyvMEiql21Wuv3w5ce47TWC7XW9VrrWq31Aq11WfinjNIR+6nbM4IN22PXaDpJGk37RUaYioTT1DKGVeu2xux4Walp7KqK3T01g4WEh0hA5bz1cYy/zB0yMVBfSXiIBORiW3Vse0AsWtPe3h7TYyY7CQ+RkGpi3OMyNr+Ij1avjukxk52Eh0hItU0mtI7NZMgAk0rLWLX8nZgdbzCQ8BAJqaE5h+qG+pgdLz8ji63yHJc+kfAQCammcQyfVMbuiW5KKRlp2kcSHiIh+QOjeTOGN8gBGANBOjslz9DUKQAADydJREFUQCIl4SESVB4fx/AGOYDRucNYu3ZtTI+ZzCQ8RIIyUB3jWcUmFZexctmymB4zmUl4iATVTnvHnpgescPvp6E+drOYJbsYTxgpRCQ0pXlX8cA1J8XsiNtrqvnDGy/zi3vujtkxk52Eh0g4uWk/557vjac4Nzcmx2tqaebWpx/jZ7//HUZjjOcSSWISHiKheJ1PcvU5QY6cMjYmx+vo7OTav/2Zm379SxwOR0yOOVhIm4dIGBbTSs458i0uPWXB/7d359FRlOkex78PaxIhBAJ4wxrkIgMioIER0eGOCCNwRMWrsigCioKODptsKgoEEAkiR0xcZkRUEAUFBGVAQBwQVBZZL8sgKIIEoyLici+gee4fVTlE7CSkk367O3k+5/TprkpVP++bTv9SVV31tpN62dnZPDTnBR5IHUtiYqKTmiWJhYeJEEfokJLOjEE9nFWc9MYcet13D8kNGjirWZJYeJgI8DOtGj/I/HF3UKaMmz/JZ/+5mMuu7cwlKSlO6pVEFh4mzLJpXHcYiyf1IC4mxknFBevXEP+HhnTs3MlJvZLKwsOEVd0aqcwadSVJidWd1Fu7cxuHyv5Kr759Cl7Y5MvCw4RNYvxLpN4ZT5uLGhe8cDHYdfAzVh7cx9+GP+CkXklnH9WasIitsIY7u+ymT6ebnNT78uss/r5uFVMz0r0raE2RWXgY50Q+o0ub2UwecJeTet//+COTFs9j6nPP2Elgxch2W4xjJ7i86VheHdPPyRbAyVOneOi1F0mdPo3Y2NiQ1ytNLDyMQ7/QNHkYb03s7eQLnbyTwGYyYsJ4qlWrFvJ6pY2Fh3Em+T/GMHfMNVRPSHBSb+L82dw++D7q1a/vpF5pY+FhnKhZ9RmeuLcezRsmO6mXsXQRV3TrSvOWLZ3UK40sPEzIVYpdzv3djnJju9ZO6s3/YDWJzZrQvmNHJ/VKKwsPE1Jly+yi258W8/DtXZ3Ue3/7Fr6KKUv33rc6qVeaWXiYEPqWdi2mMHNEbyfVdhzYz7+OfMZfhw5xUq+0s/M8TIicokXD4SxM7Uu5cqH/MzuUdZRZG9aQlj7DTgJzxMLDhIDSsNZI3hx/PVUqVQp5te9+OMHj7yxg6rMZzq7KNRYeJgSSEqfxzJCLaVi7dshr/d/Jk4x5/WUmzJhOjKOrco3HYtoUqyrnLWBkz1N0bN085LWys7N56NWZjJqUSoKjc0fMGRYepthUKPcJPa9ey6CbQv8Rqaoy/vWX6Td0EHXq1g15PfN7Fh6mmGRydUo66YN7OqmWvnQRV91yI82ah34LxwRm4WGKTOQA7VqM4o1x/ZwcsJy1chm1UlrwX+3bh7yWyZsdMDVFUrH8x3T649+ZN3agk4vdXlm9nLjGF3Bj91tCXsvkz8LDBC0+bgm9OrxHxpABTs6tmPuvVZRNrsMtt/UKeS1TMAsPE5QaCc8z5OajjL7VzWng89a+xy+1a9Dbxh6NGBYeppCUujVSmTygKr06dHFSccH6NfxUPYF+d97hpJ45NxYephBO8Z+1hzNzRCv+1KKJk4qLP/qAY1Vi6T/wbif1zLmz8DDn6AQXNxjCm+O70qhuHScV39n4IUdiyjDw3nuc1DOFYx/VmnNwmDZN7+e9J3s6C47lmzfweZlfGDjob07qmcKzLQ+Tr3LldtDh0idZkHoXsRXdXDuycssm9pz+iUEjhzupZ4Jj4WHydF7MKv673QJmjhzg7CsL3t++he3/e5yho0c5qWeCZ+FhAkqMn82ArruZ0L+vs/Ex1u7cxsfHMhn56CNO6pmiCdkxDxFpLCJbc91OiMhgEakmIitEZJ9/XzVUbTDBUGolTmVi/6NMvKubs+BYv2sHa45+wYhHxjipZ4ouZOGhqntVtaWqtgRSgJ+BhcAoYJWqNgJW+dMmIvxKg6RR/GNEAgOu+7Ozqhv27GLloU95cPxYGwUsirj6tOVqYL+qHgSuB17y578E3OCoDSZfP9Ok3r0smtCazpe5+7qCzfv28s7+XYyZOMGCI8q4OubRA5jrPz5fVTMBVDVTRGo6aoPJ09ekXDicJY/1ICmxurOqW/f/m0V7tjE+7XELjigU8i0PEakAXAfML+R6d4vIJhHZBF+HpnGGMrKP9pcMZ/X0fk6DY8eB/czbsZlxUyZbcEQpF7stnYFPVPUrf/orEUkC8O+zAq2kqs+raitVbQU1HDSz9ImpsJ4b201lWdpAKsed56zuroOfMWfrh0x4Is0GLI5iLl65npzZZQFYDORcGtkHeMtBG8xZEs5byN3XLmTe2P6Ud/DVCDn2HjrIrI1rmTjtCQuOKCeqGronF4kDDgEXqOr3/rxEYB5QD/gCuFlVj+X/PK0UNoWsnaVNzaoZjOjxHcO6X+O07qdfHubZD1Yy+anpTr7LxRSdiGz2tv5/L6SvoKr+DCSeNe9bvE9fjHPZ1Ks5jqn31OTmq9wGx4Ejh8lY8y5Tnn7KgqOEsFex1DhJozoP8PLoNrS5qLHTygePZvLU6mWkpc+w4ChB7JUsFY7T/IKhLJrQjQa1kpxWPpR1lGkrljAlfQblHYxxatyx8CjxvuCKZmNY8lgfqlaOd1r5y6+zmPLPRaRlPE3FihWd1jahZ+FRgpUvt5VrWs9g/ti7iXH85j367Tc89vabpD3ztH0NZAll4VFCVYpdTvf2S3l+2ADnH4lmfXeMCYvnMSXjaWJjY53WNu5YeJRA1avM4r5uB3i0b2/ntb85fpxxC+fyePoM4uLinNc37oT0PI/iknJpK133wYZwNyNqHDiwl6ZN3X6ikiMzM5NK8fFUrlw5LPVN8QrbeR7FRcpATJydjXiumjZzM7J5IEm1a4ettnHL3pHGmKBYeBhjgmLhYYwJioWHMSYoFh7GmKBYeBhjgmLhYYwJioWHMSYoFh7GmKBYeBhjgmLhYYwJioWHMSYoFh7GmKBYeBhjgmLhYYwJioWHMSYoFh7GmKBYeBhjghIVY5iKyA/A3nC3oxhUB74JdyOKgfUjsoSyH/VVtUagH0TFGKbA3rwGYY0mIrLJ+hE5rB9FY7stxpigWHgYY4ISLeHxfLgbUEysH5HF+lEEUXHA1BgTeaJly8MYE2EiOjxEpJOI7BWRT0VkVLjbUxgiMlNEskRkZ6551URkhYjs8++rhrONBRGRuiKyWkR2i8j/iMggf3609SNGRDaIyDa/H+P8+VHVjxwiUlZEtojI2/50WPoRseEhImWBdKAz0BToKSJNw9uqQpkFdDpr3ihglao2Alb505HsF2CYqjYB2gB/9V+DaOvHSaC9qrYAWgKdRKQN0dePHIOA3bmmw9MPVY3IG3A5sDzX9GhgdLjbVcg+JAM7c03vBZL8x0l456+EvZ2F6M9bQMdo7gcQB3wCXBaN/QDq4AVEe+Btf15Y+hGxWx5AbeBQrunD/rxodr6qZgL49zXD3J5zJiLJwCXAx0RhP/xN/a1AFrBCVaOyH8B0YASQnWteWPoRyeEhAebZR0NhICKVgDeBwap6ItztCYaq/qqqLfH+c/9RRJqFu02FJSLXAlmqujncbYHIDo/DQN1c03WAI2FqS3H5SkSSAPz7rDC3p0AiUh4vOOao6gJ/dtT1I4eqHgfexzseFW39uAK4TkQ+B14D2ovIbMLUj0gOj41AIxFpICIVgB7A4jC3qagWA338x33wjiFELBER4AVgt6pOy/WjaOtHDRFJ8B/HAh2APURZP1R1tKrWUdVkvPfDe6p6G+HqR7gPABVwcKgL8G9gP/BQuNtTyLbPBTKB03hbUXcCiXgHu/b599XC3c4C+nAl3q7idmCrf+sShf1oDmzx+7ETeMSfH1X9OKtPf+bMAdOw9MPOMDXGBCWSd1uMMRHMwsMYExQLD2NMUCw8jDFBsfAwxgTFwiOKiciTIjI41/RyEflHruknRGRoPusniMi9oW7nuRKRB4u4/g15XTwpIveLyE4RWeqfN4SIXCki0wItbwpm4RHd1gNtAUSkDN4o2hfl+nlbYF0+6ycAhQoP8YTq76ZI4QHcgHcFdiD9OXO+xzX+CXBjgNQi1iy1LDyi2zr88MALjZ3ADyJSVUQqAk2ALSJSSURWicgnIrJDRK7315kMNBSRrSKSBiAiw0Vko4hszzXuRbI/pkcG3hWpuS8bQERai8h6f7yMDSJS2R9D40W/3hYRucpftq+ILBCRZf74E1P8+ZOBWL8tc/x5t/nPt1VEnvOHaUBEfhSRiX69j0TkfBFpC1wHpPnLNwzw+yqPd1XtaaA3sFRVvyv6y1BKhftMObsV7QZ8DtQDBgAD8f6TdsG7DmKNv0w5IN5/XB34FO/Cw2R+O2TAX/DGwxS8fyxvA+385bKBNgHqVwAOAK396Xi/3jDgRX/eH4AvgBigr798FX/6IFDXX+7HXM/bBFgClPenM4Db/ccKdPUfTwEe9h/PAm7K4/fUG2+rYzZQGe9MzPLhfv2i+RYt39ti8paz9dEWmIY3bEFb4Hu83RrwwmCSiLTDC4HawPkBnusv/m2LP10JaIT3xj+oqh8FWKcxkKmqGwHUv+pWRK4EZvjz9ojIQeBCf51Vqvq9v9wuoD6/HX4B4GogBdjo7WEQy5kLvk7hBRvAZrwxRvKlqq8Ar/g1HwWeAjqLyO1+7WGqmp3PU5izWHhEv5zjHhfj7bYcwvuvfwKY6S9zK1ADSFHV0/5VmTEBnkuAx1T1ud/M9Mby+CmP+kLgoRICDamQ42Sux78S+O9QgJdUdXSAn51Wf3Min/UDEpFaeFtJ40RkA96gUxPxwmrFuT6PsWMeJcE64FrgmHpjVhzDOxB6OfChv0wVvHEgTvvHHur783/A24TPsRy4wx+/AxGpLSIFDSyzB6glIq39dSqLSDlgDV5oISIX4u1aFfSVoaf9IQDA2624Kae+P05n/bxXDdifQFLxDpSCtzWjeFtjcQWsZ85i4RH9duAdx/jorHnfq2rO95fOAVqJyCa8N/QeAFX9Fljnf4SZpqrvAq8CH4rIDuANCngzquopoDswQ0S24f33jsE7RlHWf57Xgb6qejLvZwK84y3bRWSOqu4CHgbeFZHt/vMmFbD+a8Bw/wDt7w6YisglfptzdstewPtdXQosK+C5zVnsqlpjTFBsy8MYExQLD2NMUCw8jDFBsfAwxgTFwsMYExQLD2NMUCw8jDFBsfAwxgTl/wGiwGvPVFrr3wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 288x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# PLot an SNMR water content profile\n",
    "# Finally lets plot an inversion\n",
    "\n",
    "fig, ax1 = plt.subplots(1,1, sharey=True, figsize = (4,8))\n",
    "\n",
    "inv_id = 40\n",
    "\n",
    "ax1 = SNMR_utils.plot_profile(ax1, df_inversions[df_inversions['inversion_id']==inv_id],\n",
    "                              doi = None,  plot_mobile_water = True)\n",
    "\n",
    "plt.gca().invert_yaxis()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   borehole_id     log_K     Easting     Northing Borehole_name  \\\n",
      "0       635728  1.711760  520144.760  8306293.880      RN040873   \n",
      "1       635735 -0.102018  505309.448  8293856.701      RN040878   \n",
      "2       635737  0.134896  507769.229  8299897.749      RN040884   \n",
      "3       635740  2.417096  500053.109  8296644.456      RN040882   \n",
      "4       635741  1.860120  500049.774  8295302.212      RN040876   \n",
      "5       635742  1.206960  500054.317  8293767.687      RN040875   \n",
      "\n",
      "  Alternative_name  Ground_elevation_mAHD  Depth_from  Depth_to  \\\n",
      "0             KR08                  7.519          31        37   \n",
      "1             KR33                 16.494          36        42   \n",
      "2             KR38                 12.462          29        35   \n",
      "3             KR46                 18.357          17        23   \n",
      "4             KR48                 19.046          20        26   \n",
      "5             KR49                 19.951          29        35   \n",
      "\n",
      "  Construction_name         stratigraphy  \n",
      "0            screen   Cenozoic sediments  \n",
      "1            screen  Permo-Carboniferous  \n",
      "2            screen  Permo-Carboniferous  \n",
      "3            screen   Cenozoic sediments  \n",
      "4            screen   Cenozoic sediments  \n",
      "5            screen   Cenozoic sediments  \n"
     ]
    }
   ],
   "source": [
    "# Bring in the slug test processed K values\n",
    "\n",
    "infile = r\"C:\\Users\\PCUser\\Desktop\\EK_data\\Boreholes\\EK_slug_test_K_processed.csv\"\n",
    "\n",
    "df_slug = pd.read_csv(infile)\n",
    "\n",
    "print(df_slug)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we do a spatial query to find the representative SNMR observations\n",
    "\n",
    "\n",
    "SNMR_coords = df_acquisitions[['mid_X', 'mid_Y']].values\n",
    "\n",
    "distances, indices = spatial_functions.nearest_neighbours(df_slug[['Easting', 'Northing']],\n",
    "                                                         SNMR_coords,\n",
    "                                                         points_required = 1,\n",
    "                                                         max_distance = 200.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_coloc_bores = df_slug.copy()\n",
    "\n",
    "df_coloc_bores['colocated_SNMR_acqu_id'] = -999\n",
    "\n",
    "# Iterate through the bores and if a site was returned\n",
    "# add the site id to the dataframe\n",
    "for i, (index, row)  in enumerate(df_coloc_bores.iterrows()):\n",
    "    idx = indices[i]\n",
    "    if idx != len(df_acquisitions):\n",
    "        # Since we only have one acquisition per site we can use site_id\n",
    "        df_coloc_bores.at[index, 'colocated_SNMR_acqu_id'] = idx\n",
    "        \n",
    "# Remove the bores with no location\n",
    "\n",
    "df_coloc_bores = df_coloc_bores[df_coloc_bores['colocated_SNMR_acqu_id'] != -999]\n",
    "\n",
    "# Get the SNMR site_ids\n",
    "\n",
    "SNMR_ids = df_coloc_bores['colocated_SNMR_acqu_id'].unique()\n",
    "\n",
    "# Trim the inversions\n",
    "\n",
    "df_inversions_subset = df_inversions[df_inversions['acquisition_id'].isin(SNMR_ids)]\n",
    "\n",
    "df_inversions_subset=  df_inversions_subset.merge(df_acquisitions[df_acquisitions.index.isin(SNMR_ids)], \n",
    "                                              how='inner', left_on='acquisition_id',\n",
    "                                              right_index=True, )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>borehole_id</th>\n",
       "      <th>log_K</th>\n",
       "      <th>Easting</th>\n",
       "      <th>Northing</th>\n",
       "      <th>Borehole_name</th>\n",
       "      <th>Alternative_name</th>\n",
       "      <th>Ground_elevation_mAHD</th>\n",
       "      <th>Depth_from</th>\n",
       "      <th>Depth_to</th>\n",
       "      <th>Construction_name</th>\n",
       "      <th>stratigraphy</th>\n",
       "      <th>colocated_SNMR_acqu_id</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>635728</td>\n",
       "      <td>1.711760</td>\n",
       "      <td>520144.760</td>\n",
       "      <td>8306293.880</td>\n",
       "      <td>RN040873</td>\n",
       "      <td>KR08</td>\n",
       "      <td>7.519</td>\n",
       "      <td>31</td>\n",
       "      <td>37</td>\n",
       "      <td>screen</td>\n",
       "      <td>Cenozoic sediments</td>\n",
       "      <td>223</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>635735</td>\n",
       "      <td>-0.102018</td>\n",
       "      <td>505309.448</td>\n",
       "      <td>8293856.701</td>\n",
       "      <td>RN040878</td>\n",
       "      <td>KR33</td>\n",
       "      <td>16.494</td>\n",
       "      <td>36</td>\n",
       "      <td>42</td>\n",
       "      <td>screen</td>\n",
       "      <td>Permo-Carboniferous</td>\n",
       "      <td>213</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>635737</td>\n",
       "      <td>0.134896</td>\n",
       "      <td>507769.229</td>\n",
       "      <td>8299897.749</td>\n",
       "      <td>RN040884</td>\n",
       "      <td>KR38</td>\n",
       "      <td>12.462</td>\n",
       "      <td>29</td>\n",
       "      <td>35</td>\n",
       "      <td>screen</td>\n",
       "      <td>Permo-Carboniferous</td>\n",
       "      <td>217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>635740</td>\n",
       "      <td>2.417096</td>\n",
       "      <td>500053.109</td>\n",
       "      <td>8296644.456</td>\n",
       "      <td>RN040882</td>\n",
       "      <td>KR46</td>\n",
       "      <td>18.357</td>\n",
       "      <td>17</td>\n",
       "      <td>23</td>\n",
       "      <td>screen</td>\n",
       "      <td>Cenozoic sediments</td>\n",
       "      <td>215</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>635741</td>\n",
       "      <td>1.860120</td>\n",
       "      <td>500049.774</td>\n",
       "      <td>8295302.212</td>\n",
       "      <td>RN040876</td>\n",
       "      <td>KR48</td>\n",
       "      <td>19.046</td>\n",
       "      <td>20</td>\n",
       "      <td>26</td>\n",
       "      <td>screen</td>\n",
       "      <td>Cenozoic sediments</td>\n",
       "      <td>127</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>635742</td>\n",
       "      <td>1.206960</td>\n",
       "      <td>500054.317</td>\n",
       "      <td>8293767.687</td>\n",
       "      <td>RN040875</td>\n",
       "      <td>KR49</td>\n",
       "      <td>19.951</td>\n",
       "      <td>29</td>\n",
       "      <td>35</td>\n",
       "      <td>screen</td>\n",
       "      <td>Cenozoic sediments</td>\n",
       "      <td>253</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   borehole_id     log_K     Easting     Northing Borehole_name  \\\n",
       "0       635728  1.711760  520144.760  8306293.880      RN040873   \n",
       "1       635735 -0.102018  505309.448  8293856.701      RN040878   \n",
       "2       635737  0.134896  507769.229  8299897.749      RN040884   \n",
       "3       635740  2.417096  500053.109  8296644.456      RN040882   \n",
       "4       635741  1.860120  500049.774  8295302.212      RN040876   \n",
       "5       635742  1.206960  500054.317  8293767.687      RN040875   \n",
       "\n",
       "  Alternative_name  Ground_elevation_mAHD  Depth_from  Depth_to  \\\n",
       "0             KR08                  7.519          31        37   \n",
       "1             KR33                 16.494          36        42   \n",
       "2             KR38                 12.462          29        35   \n",
       "3             KR46                 18.357          17        23   \n",
       "4             KR48                 19.046          20        26   \n",
       "5             KR49                 19.951          29        35   \n",
       "\n",
       "  Construction_name         stratigraphy  colocated_SNMR_acqu_id  \n",
       "0            screen   Cenozoic sediments                     223  \n",
       "1            screen  Permo-Carboniferous                     213  \n",
       "2            screen  Permo-Carboniferous                     217  \n",
       "3            screen   Cenozoic sediments                     215  \n",
       "4            screen   Cenozoic sediments                     127  \n",
       "5            screen   Cenozoic sediments                     253  "
      ]
     },
     "execution_count": 170,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_coloc_bores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For each bore with slug test results we will find all of the SNMR observations\n",
    "# for the screened interval and average them and add them to the merge dataframe\n",
    "\n",
    "# Define the columns to add to the dataframe\n",
    "cols = ['Mobile_water_content', 'Bound_water_content',\n",
    "        'Total_water_content', 'T2*']\n",
    "\n",
    "for c in cols:\n",
    "    df_coloc_bores[c] = np.nan\n",
    "\n",
    "for index, row in df_coloc_bores.iterrows():\n",
    "    \n",
    "    # Subset the interpreted dataframe\n",
    "    SNMR_key = row.colocated_SNMR_acqu_id\n",
    "    \n",
    "   \n",
    "    # Get the stratigraphy for the borehole and average\n",
    "    df_temp = df_inversions_subset[df_inversions_subset['acquisition_id'] == SNMR_key]\n",
    "    \n",
    "    # Subset based on the screenen interval and average\n",
    "    \n",
    "    mask = (df_temp['Depth_from'] >= row.Depth_from) & \\\n",
    "           (df_temp['Depth_to'] <= row.Depth_to)\n",
    "    \n",
    "    interval_SNMR = df_temp[mask].mean()\n",
    "    \n",
    "    # Add this to the df_merged dataframe\n",
    "    df_coloc_bores.at[index, cols] = interval_SNMR[cols]\n",
    "    \n",
    "# Remove any nulls where the screen didn't intersect the layered model\n",
    "df_coloc_bores.dropna(how=\"any\", subset = cols, inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEGCAYAAAC3lehYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3yV5f3/8dcnCQl7bwh7g4gSAcE9cNdZxUkVpVX7U6u2lbbfWkfr+HZY9atV3FRRHFWK4ihqxQEIspeEHVaAMBLCyPj8/jh39BghHElOzsj7+Xicx7nPdc59n88FIR+ucV+XuTsiIiLRkBLrAEREJHkpyYiISNQoyYiISNQoyYiISNQoyYiISNSkxTqAeNG8eXPv1KlTrMMQEUkos2bN2uLuLQ70vpJMoFOnTsycOTPWYYiIJBQzW13R++ouExGRqFGSERGRqFGSERGRqFGSERGRqFGSERGRqFGSERGRqFGSERGRqNF9MiIiNdibs9fhOOcNaIeZVfn1o9aSMbNnzCzXzBaElf3YzBaaWamZZZX7/BgzyzazpWZ2Wlj5QDObH7z3sAV/CmaWYWavBOXTzaxT2DkjzWxZ8BgZrTqKiCSy/D1F3DNpEa98uTZq3xHN7rLngNPLlS0ALgA+CS80sz7ACKBvcM5jZpYavP04MBroHjzKrjkK2Obu3YC/AQ8E12oK3AkMBgYBd5pZk6qsmIhIMhj7yQq27trHmDN6R6UVA1FMMu7+CZBXrmyxuy/dz8fPBV52973uvhLIBgaZWRugobt/4aEtPF8Azgs75/ng+DXg5KCVcxrwgbvnufs24AO+n+xERGq03J17GDt1JWf1b8PhmY2j9j3xMvDfDghvr+UEZe2C4/Ll3znH3YuBHUCzCq71PWY22sxmmtnMzZs3V0E1REQSw0NTllFUUsovh/eM6vfES5LZXzvNKyg/1HO+W+j+pLtnuXtWixYHXERURCSpLN9cwCtfruXywR3o1LxeVL8rXpJMDpAZ9ro9sD4ob7+f8u+cY2ZpQCNC3XMHupaIiAAPvruE2mkp/L+Tu0f9u+IlyUwERgQzxjoTGuCf4e4bgHwzGxKMt1wFvBV2TtnMsYuAD4Nxm/eA4WbWJBjwHx6UiYjUeLNW5/Hewk389PiuNK+fEfXvi9p9MmY2HjgBaG5mOYRmfOUBjwAtgLfNbI67n+buC81sArAIKAZudPeS4FLXE5qpVgeYHDwAngbGmVl2cN0RAO6eZ2b3AF8Gn7vb3b8zAUFEpCZyd+57ZwktGmRw7bGdq+U7LfSff8nKynJtWiYiyez9hRsZPW4Wfzy/H5cP7lgl1zSzWe6edaD346W7TEREoqiopJT7Jy+ha4t6XJKVefATqoiSjIhIDfDitNWs2LKL35zZm7TU6vvVryQjIpLkduwu4u9TljG0azNO6tWyWr9bSUZEJMn930fZbN9dxG/Pit7yMQeiJCMiksTW5hXy3GeruPDI9vRt26jav19JRkQkid3/7hJSU4zbo7x8zIEoyYiIJKlZq7fx9rwNXHdcF1o3qh2TGJRkRESSkLtz79uLaNEgg58e1yVmcSjJiIgkoUnzNjB7zXZ+Obwn9TJitwmykoyISJLZU1TCA+8uoVfrBlw4sP3BT4giJRkRkSTz/OeryNm2m9+d1YfUlOqdslyekoyISBLJ27WPRz/K5sSeLTime/NYh6MkIyKSTB76z9cU7ivhN2f2jnUogJKMiEjSWLoxnxenr+HSQZl0b9Ug1uEASjIiIknB3bl70kLqZ6Rx26mxufFyf5RkRESSwPuLNvFZ9lZ+cUp3mtRLj3U431CSERFJcHuKSvjj24vp0ao+Vwypms3IqoqSjIhIgnvms5WsySvk92f3rda9YiIRX9GIiMgPsmnnHh79MJtT+7SKiynL5SnJiIgksAfeXUJxifO7s+JjynJ5SjIiIglq9pptvPHVOkYd25mOzerFOpz9UpIREUlApaXOH/69iJYNMrjxxG6xDueAlGRERBLQv2avY+7a7fz69F7Uj+EqywejJCMikmAK9hbzwLtLODyzMecf0S7W4VRISUZEJMH830fZ5Obv5Q/n9CElxqssH4ySjIhIAlmxuYCnpq7ggiPbcUSHJrEO56CUZEREEoS7c+fEhdROS2XMGfE5Zbk8JRkRkQTx7oKNTF22hVuH96BFg4xYhxMRJRkRkQRQuK+YeyYtolfrBlwZZ+uTVeSgScbMbjazhhbytJl9ZWbDqyM4EREJefTDbNbv2MM95/WLu/XJKhJJpNe4+05gONACuBq4P6pRiYjIN5ZvLmDs1BVceGR7jurUNNbh/CCRJJmy+XFnAs+6+9ywMhERiSJ35w8TF1K7Vip3nNEr1uH8YJEkmVlm9j6hJPOemTUASqMbloiIwLeD/bedmjiD/eEiWYtgFDAAWOHuhWbWjFCXmYiIRFHhvmLunrSI3m0axt1mZJE6YJIxsyPLFXUxUy+ZiEh1eeTDbDbs2MMjlx6RUIP94SpqyfylgvccOKmKYxERkcDy4M7+iwa2JyvBBvvDHTDJuPuJ1RmIiIiEuDu/f2tBwg72h4tofWgz6wf0AWqXlbn7C9EKSkSkJntrzno+y97KPef2pXn9xBvsDxfJzZh3Ao8EjxOBB4EfRXDeM2aWa2YLwsqamtkHZrYseG4S9t4YM8s2s6VmdlpY+UAzmx+897AFA0NmlmFmrwTl082sU9g5I4PvWGZmIyP6kxARiQPbC/dxz6RFDMhszOWDE3OwP1wkI0kXAScDG939auBwIJLU+hxwermyO4Ap7t4dmBK8xsz6ACOAvsE5j5lZanDO48BooHvwKLvmKGCbu3cD/gY8EFyrKXAnMBgYBNwZnsxEROLZ/ZOXsH13EX86/7C4X8Y/EpEkmd3uXgoUm1lDIBfocrCT3P0TIK9c8bnA88Hx88B5YeUvu/ted18JZAODzKwN0NDdv3B3B14od07ZtV4DTg5aOacBH7h7nrtvAz7g+8lORCTufLkqj5e/XMu1x3SmT9uGsQ6nSkQyJjPTzBoDY4FZQAEw4xC/r5W7bwBw9w1m1jIobwdMC/tcTlBWFByXLy87Z21wrWIz2wE0Cy/fzznfYWajCbWS6NChwyFWSUSk8vYVl/KbN+bTrnEdbj6le6zDqTIHTTLufkNw+A8ze5dQy2JeFcexvzahV1B+qOd8t9D9SeBJgKysrP1+RkSkOoyduoJluQU8PTKLuukRzclKCJEM/B9X9gA6AI2D40OxKegCI3jODcpzgMywz7UH1gfl7fdT/p1zzCwNaESoe+5A1xIRiUurt+7i4SnLOPOw1pzcu1Wsw6lSkYzJ/DLs8T/Av4E/HOL3TQTKZnuNBN4KKx8RzBjrTGiAf0bQtZZvZkOC8Zaryp1Tdq2LgA+DcZv3gOFm1iQY8B8elImIxB1353dvLqBWagp3ntM31uFUuUi6y84Jf21mmYSmMVfIzMYDJwDNzSyH0Iyv+4EJZjYKWAP8OPiOhWY2AVgEFAM3untJcKnrCc1UqwNMDh4ATwPjzCybUAtmRHCtPDO7B/gy+Nzd7l5+AoKISFyYOHc9U5dt4a4f9aVVw9oHPyHBWOg//z/ghFCLYp67HxadkGIjKyvLZ86cGeswRKQG2VFYxMl//S/tGtfmjRuGkZqAU5bNbJa7Zx3o/YO2ZMzsEb4dOE8htCLz3KoJT0Sk5rr/3cXk7drLc1cflZAJJhIRTWEOOy4Gxrv7Z1GKR0SkRvh8+RbGz1jL6OO60K9do1iHEzWRjMk8f7DPiIhI5HbvK2HMG/Pp2KwuvzilR6zDiaqK9pOZzwHuLwFw9/5RiUhEJMk99J+vWb21kJeuG0yd9NSDn5DAKmrJnB083xg8jwueLwcKoxaRiEgSm5eznbFTV3DpoEyGdm0e63CirqL9ZFYDmNkwdx8W9tYdZvYZcHe0gxMRSSb7ikv51WvzaNEggzvO6B3rcKpFJDdj1jOzY8pemNlQoF70QhIRSU5P/Hc5Szbmc+95h9GoTq1Yh1MtIpldNgp4xszKpj9sB66JXkgiIsknOzefRz7M5uz+bTi1T3ItHVORSGaXzQIOD5b5N3ffEf2wRESSR0mp86vX5lE3I5U//Cj5lo6pSEWzy65w93+a2a3lygFw979GOTYRkaQw7otVfLVmO3+9+PCE3075h6qoJVM27tKgOgIREUlGa/MKefC9pRzfowXnH7Hfra2SWkWzy54Inu+qvnBERJJHaanz69fnYcAfz+/3TU9QTRLJfjIPmllDM6tlZlPMbIuZXVEdwYmIJLIXp6/m8+Vb+e1ZfWjfpG6sw4mJSKYwD3f3nYRuzswBehDaW0ZERA5gzdZC/vTOEo7t3pxLB2Ue/IQkFUmSKZvMfSahxTG1N4uISAVKS53bX5tLWorxwIX9a2Q3WZlI7pP5t5ktAXYDN5hZC2BPdMMSEUlcz32+ihkr8/jfi/rTtnGdWIcTUwdtybj7HcDRQJa7FxFat+zcaAcmIpKIVmwu4MH3lnBSr5ZcNLB9rMOJuUgG/usSWiTz8aCoLXDAXdBERGqqklLn9lfnkpGWyn0XHFaju8nKRDIm8yywDxgavM4B7o1aRCIiCerpT1fw1Zrt3PWjvrRqWDvW4cSFSJJMV3d/ECgCcPfdgNKziEiYZZvy+fP7XzO8TyvOHdA21uHEjUiSzD4zq0OwgZmZdQX2RjUqEZEEUlxSyu2vzqVeeip/PF/dZOEimV12J/AukGlmLwLDgJ9EMygRkUTy6EfZzM3ZwaOXHUGLBjVrbbKDqTDJWCgdLwEuAIYQ6ia72d23VENsIiJxb/aabTzyYTbnH9GOs/urm6y8CpOMu7uZvenuA4G3qykmEZGEsGtvMb94ZQ6tG9bmrnNr1hL+kYpkTGaamR0V9UhERBLMvW8vZnVeIX+9+HAa1q4ZO13+UJGMyZwI/NTMVgO7CHWZubv3j2pkIiJx7INFmxg/Yw0/Pb4Lg7s0i3U4cSuSJHNG1KMQEUkgm/P3csfr8+jTpiG3ntoj1uHEtUi2X15dHYGIiCQCd+eO1+eRv7eY8SMGkJGWGuuQ4lokYzIiIhJ4acYapizJZcwZvejRShsHH4ySjIhIhFZsLuDeSYs5tntzRh7dKdbhJAQlGRGRCOwrLuWWV+aQUSuFP//4cFJSdFd/JA44JmNm+QRLyeyPuzeMSkQiInHof99bwrycHfzjioFa/PIHOGCScfcGAGZ2N7ARGEdo+vLlgDoiRaTG+GhpLmOnruTKIR05vV/rWIeTUCLpLjvN3R9z93x33+nujwMXRjswEZF4kLtzD7dPmEuv1g347Vm9Yx1OwokkyZSY2eVmlmpmKWZ2OVAS7cBERGKttNS5dcJcdu0r5pFLj6B2LU1X/qEiSTKXARcDm4LHj4MyEZGk9o9PlvNp9hbuPKcv3TVd+ZBEcjPmKuDc6IciIhI/vlqzjb+8/zVnHdaGEUdlxjqchHXQloyZ9TCzKWa2IHjd38x+F/3QRERiY8fuIm4aP5vWDWvzpwu0CVllRNJdNhYYw7fbL88DRlTmS83sZjNbYGYLzeyWoKypmX1gZsuC5yZhnx9jZtlmttTMTgsrH2hm84P3Hg72v8HMMszslaB8upl1qky8IlJzuDu/+dd8NuzYw8OXHkGjOlpduTIiSTJ13X1GubLiQ/1CM+sHXAcMAg4Hzjaz7sAdwBR37w5MCV5jZn0IJbW+wOnAY2ZWNvr2ODAa6B48Tg/KRwHb3L0b8DfggUONV0Rqlhenr+HteRu49dQeDOzY5OAnSIUiSTJbzKwrwY2ZZnYRsKES39kbmObuhe5eDPwXOJ/QuM/zwWeeB84Ljs8FXnb3ve6+EsgGBplZG6Chu3/h7g68UO6csmu9Bpxc1soRETmQ+Tk7uPvfiziuRwuuP75rrMNJCpEkmRuBJ4BeZrYOuAX4WSW+cwFwnJk1M7O6wJlAJtDK3TcABM8tg8+3A9aGnZ8TlLULjsuXf+ecIJHtAL634YOZjTazmWY2c/PmzZWokogkuh2FRdzw0iya1U/noUsGaNmYKlLh7DIzSwGy3P0UM6sHpLh7fmW+0N0Xm9kDwAdAATCXirvf9vc37RWUV3RO+VieBJ4EyMrKOuASOiKS3Nyd21+by4bte3jlp0fTtF56rENKGhW2ZNy9FPh5cLyrsgkm7LpPu/uR7n4ckAcsAzYFXWAEz7nBx3MItXTKtAfWB+Xt91P+nXPMLA1oFHyPiMj3jJ26gg8WbWLMmb01DlPFIuku+8DMbjezzGAGWFMza1qZLzWzlsFzB+ACYDwwERgZfGQk8FZwPBEYEcwY60xogH9G0KWWb2ZDgvGWq8qdU3ati4APg3EbEZHv+HJVHg+8u5Qz+rXmmmGdYh1O0olk++Vrgucbw8oc6FKJ733dzJoRmhZ9o7tvM7P7gQlmNgpYQ2hlAdx9oZlNABYR6la70d3LlrW5HngOqANMDh4ATwPjzCybUAumUlOuRSQ5bSnYy89f+orMJnV44KL+uh8mCkz/wQ/JysrymTNnxjoMEakmJaXOyGdmMGNVHv+6YSh92zaKdUgJycxmuXvWgd6PpCVTdm9LH+CbTRTc/YXKhyciEht/n7KMT7O3cP8FhynBRNFBk4yZ3QmcQCjJvAOcAXxK6L4UEZGE88GiTTw8ZRkXDWzPJVqXLKoiGfi/CDgZ2OjuVxO6Sz8jqlGJiETJ8s0F/OKVOfRv34h7z+uncZgoiyTJ7A6mMhebWUNCU4srM+gvIhIT+XuKGP3CTDLSUvjHFQO1P0w1iGRMZqaZNSa0UOYsQjdQll/LTEQkrpWWOrdNmMuqrYX8c9Rg2jauE+uQaoRI9pO5ITj8h5m9S2i9sHnRDUtEpGo99nE27y/axO/P7sPRXb+3ypRESSQD/8ftr8zdP4lOSCIiVeujpbn85YOvOW9AW67WDZfVKpLusl+GHdcmtET/LOCkqEQkIlKFVm3Zxc3jZ9O7dUPuu0A3XFa3SLrLzgl/bWaZwINRi0hEpIoU7C3mp+NmkZJiPHHlQOqka6C/ukV0M2Y5OUC/qg5ERKQqlZQ6N4+fTfbmAl64ZhCZTevGOqQaKZIxmUf4dpn8FGAAoeX5RUTi1oPvLmHKklzuObcvw7o1j3U4NVZEU5jDjouB8e7+WZTiERGptFdnruWJT1Zw5ZCOXHl0p1iHU6NFMibz/ME+IyISL2auyuO3/1rAsG7N+P05fWIdTo0XSXfZfPazqySh3Sfd3ftXeVQiIodgbV4hPx03i3ZN6vDYZQOplRrJoiYSTZF0l5Xt0TIueL4cKATUwhGRuFGwt5jrXphJUUkpT43MolHdWrEOSYgsyQxz92Fhr+8ws8/c/e5oBSUi8kOUlDq3vDyHZbkFPHf1UXRtUT/WIUkgkrZkPTM7puyFmQ0F6kUvJBGRH+aPby/mP4tDS8Yc271FrMORMJG0ZEYBz5hZ2a4+2/l2S2YRkZh69rOVPPPZSq4e1omRQzvFOhwpJ5LZZbOAw4Nl/s3dd0Q/LBGRg3tv4UbunrSI0/q24ndnaSZZPDpod5mZ3RwkmHzgL2b2lZkNj35oIiIHNnvNNm4aP5vD2zfmoUuOIDVFa5LFo0jGZK5x953AcKAlcDVwf1SjEhGpwOqtu7j2+Zm0alibp0ZmaU2yOBZJkin778GZwLPuPjesTESkWm3btY+rn/2SEneeu/oomtfXbvDxLJIkM8vM3ieUZN4zswZAaXTDEhH5vj1FJYweN5Oc7bsZe1UWXTRVOe5FOrtsALDC3QvNrBmhLjMRkWpTdi/Ml6u28ehlR3BUp6axDkkiEMnsslLgq7DXW4Gt0QxKRCScu/O7N+fz7sKN/P7sPpzdv22sQ5IIaWEfEYl7f3n/a8bPWMuNJ3blmmM6xzoc+QEOmGTMTH+TIhJzz3y6kkc/yubSQZncPrxnrMORH6iilsxrAGY2pZpiERH5jjdnr+PuSYs4vW9r7j3vMMw0sTXRVDQmk2JmdwI9zOzW8m+6+1+jF5aI1HQfL83l9lfnMqRLUx4aMUA3WyaoiloyI4A9hBJRg/08RESiYtbqbVz/z6/o2boBY6/KonYt3WyZqA7YknH3pcADZjbP3Scf6HMiIlVpwbod/OTZGbRqmMFzVw+iQW3tC5PIIpld9rmZ/dXMZgaPv4StyCwiUmWWbszniqen07B2LV68bggtGuhu/kQXSZJ5htDimBcHj53As9EMSkRqnuWbC7j8qelkpKXw0nWDade4TqxDkioQyR3/Xd39wrDXd5nZnGgFJCI1z9q8Qi4fOx1358Vrj6ZjM+2LmCwiacnsLrcz5jBgd/RCEpGaZMOO3Vw6dhp7ikv457WD6dZS65Elk0haMj8DXggbh9kGjIxeSCJSU+Tm7+HysdPZUVjEi9cNpnebhrEOSapYJGuXzeXbnTEJ9pYREamU3J17uHTsNDbu3MML1wyif/vGsQ5JoiCSlgyg5CIiVWfjjj1cNnYam3bu4bmrB5GlFZWTVkwWyDSzX5jZQjNbYGbjzay2mTU1sw/MbFnw3CTs82PMLNvMlprZaWHlA81sfvDewxasOWFmGWb2SlA+3cw6VX8tRWR/NuzYzYgnvyA3fy/PXzOIQZ2VYJJZtScZM2sH3ARkuXs/IJXQ6gJ3AFPcvTswJXiNmfUJ3u8LnA48ZmZlt/8+DowGugeP04PyUcA2d+8G/A14oBqqJiIHsW77bi55YhpbC/bx/DVqwdQEB+0uC36hnwV0Cv98JdcuSwPqmFkRUBdYD4wBTgjefx74GPg1cC7wsrvvBVaaWTYwyMxWAQ3d/YsgzheA84DJwTl/CK71GvComZm7eyViFpFKyNlWyKVjp7G9sIhx1w5mQKbGYGqCSMZk/k1oDbP5VMG2y+6+zsz+DKwhNBX6fXd/38xaufuG4DMbzKxlcEo7YFrYJXKCsqLguHx52Tlrg2sVm9kOoBmwJTwWMxtNqCVEhw4dKls1ETmANVtDCSZ/TxEvXjtYg/w1SCRJpr2796+qLwzGWs4FOgPbgVfN7IqKTtlPmVdQXtE53y1wfxJ4EiArK0utHJEo+HpTPlc8NZ19JaW8dN0Q+rXTqlQ1SSRjMpPNbHgVfucpwEp33+zuRcAbwFBgk5m1AQiec4PP5wCZYee3J9S9lhMcly//zjlmlgY0AvKqsA4iEoG5a7dz8RNfAPDK6KOVYGqgSJLMNOBfZrbbzHaaWb6ZVWY68xpgiJnVDWaDnQwsBiby7U2eI4G3guOJwIhgxlhnQgP8M4KutXwzGxJc56py55Rd6yLgQ43HiFSvL5Zv5bKx02hQO43XfjaUnq21Q0hNFEl32V+Ao4H5VfGL2t2nm9lrwFdAMTCbUJdVfWCCmY0ilIh+HHx+oZlNABYFn7/R3UuCy10PPAfUITTgX7YlwdPAuGCSQB6h2WkiUk3+s2gTN7z0FR2b1uWf1w6mVcPasQ5JYsQOljfM7D3gDHev9KB/PMvKyvKZM2fGOgyRhPfWnHXcOmEu/do25LmrB9GkXnqsQ5IoMrNZ7p51oPcjaclsAD42s8nA3rJCbb8sIuU98+lK7nl7EUM6N2PsyCzqZ0S8qIgkqUh+AlYGj/TgISLyHaWlzp/eWcxTn67ktL6t+PuII7RlsgCRLZB5V3UEIiKJaW9xCbdNmMukeRv4ydBO/M/ZfUhN2d9dBFITRXLH/0fs/x6Tk6ISkYgkjB2FRYweN5PpK/MYc0YvRh/XhWAJQREgsu6y28OOawMXEprlJSI12Lrtu/nJMzNYtXUXfx8xgHMHtDv4SVLjRNJdNqtc0Wdm9t8oxSMiCWBeznaue2EmhXtLeP6aQQzt2jzWIUmciqS7LHyZ1BRgINA6ahGJSFx7e94Gbnt1Ds3qZfDq9YPo1Vq7WcqBRdJdNotv1worJjTTbFQ0gxKR+OPuPDwlm7/952sGdmzCE1cOpHn9jFiHJXEuku6yztURiIjErz1FJfzytXn8e+56LjiyHfddcBgZaZqiLAd3wLXLzOwoM2sd9voqM3sr2IFSOw2Fyc3fE+sQRKImd+ceLnniCybNW88dZ/TiLz8+XAlGIlbRAplPAPsAzOw44H7gBWAHwfL4Eton4+Q//5c/vr2I4pKkXnlHaqBZq/M4+5FPWZZbwJNXZvGz47tqirL8IBUlmVR3L1se/xLgSXd/3d3/B+gW/dASQ+tGtbngyHaMnbqSy5+azub8vQc/SSTOuTvPf76KS56YRp30VN64YSin9mkV67AkAVWYZIK9WCC0HP+HYe9pQaJAeloKd53bj79dcjhzc7Zz9iNTmbV6W6zDEjlku/eVcOuEudw5cSHH92jBxJ8foxlkcsgqSjLjgf+a2VuEtkmeCmBm3Qh1mUmY849ozxvXDyMjLZURT37BuC9WoS1sJNGs3rqL8x/7jDfnrOO2U3sw9qosGtWpFeuwJIFVuNS/mQ0B2gDvu/uuoKwHUN/dv6qeEKtHVS31v6OwiF9MmMOHS3K54Ih23Ht+P+qmq+En8e+DRZu4bcIczIy/jxjACT1bxjokSQCVWurf3aftp+zrqggsWTWqW4unrsrikQ+zeWjK18zJ2c6jlx5Jn7bqbpD4tLe4hPveWcJzn6+iX7uGPH75QDKb1o11WJIkItl+WX6glBTj5lO68+KowRTsKea8xz7j+c/VfSbxZ8XmAi547HOe+3wV1wzrzOvXD1WCkSqlJBNFQ7s1Z/LNxzKsazPunLiQ0eNmsb1wX6zDEgHgja9yOPuRT1m3fTdPXZXF78/po/tfpMopyURZs/oZPD3yKH53Vm8+XprLGX+fyvQVW2MdltRg+XuKuHXCnGCL5EZMvvlYTtH0ZIkSJZlqkJJiXHtsl2D2WQojxk7jT+8sZk9RSaxDkxrm8+VbOP2hqbw5ex03ndydl64bTJtGdWIdliQxJZlqdFj7Rrx907FcOqgDT36ygnMe+ZT5OZoNLtG3p6iEu/+9iMvGTic9LYXXrh/Kraf2IC1VvwIkuvQTVs3qZaTxp/MP47mrj2LnniLOf+wz/v6fZRRpSRqJknk52znr4ak889lKRh7dkbdvOoYjOzSJdVhSQyjJxMgJPVvy/i3Hc1b/NvztP15sGvsAAA7GSURBVF9z4eOfs3RjfqzDkiSyp6iEP7+3lPMf+5xde0sYN2oQd52r+7akelV4M2ZNUlU3Yx6Kd+Zv4HdvLmDn7iJ+dnxXfn5SN2rX0iwfOXTTVmzlN2/MZ8WWXVx4ZHt+f04f3bkvUVGpmzGlepx5WBuGdGnGvW8v4tGPsnln/gb+dMFhDOnSLNahSYLZsbuI+ycvYfyMNWQ2rcO4UYM4tnuLWIclNZhaMoFYtmTCfbpsC7/513zW5BVySVYmvzmzN43q6n+gUjF3590FG7lz4kK2FOzl2mO7cMsp3dU1JlF3sJaMkkwgXpIMhFbB/fuUZYyduoLGdWrxq9N78uOBmaSkaB8P+b7s3ALu+vdCpi7bQp82DXngwv4c1r5RrMOSGkJJJkLxlGTKLFq/kzsnLuDLVdvo374Rf/hRX80Kkm/k7yni4SnLePazVdRJT+W2U3twxZCOmpYs1UpJJkLxmGQg1A0yce56/vTOYjbt3MuFR7bn16f3pGXD2rEOTWKktNR5c8467pu8hC0Fe7l4YCa/PL0nzetnxDo0qYE08J/gzIxzB7TjlN6tePSjbJ6eupJ3F2zgZ8d3ZdSxndXnXsN8nr2F+yYvYf66HRzevhFjr8piQGbjWIclckBqyQTitSVT3sotu7jvncW8v2gTLRpkcMsp3bk4K5Na6iJJaks27uT+yUv4eOlm2jaqzW3De3L+Ee00Ticxp+6yCCVKkikza3Ue972zhJmrt9GleT1+eVpPTu/XGjP90kkm67bv5m8ffM3rX+XQICONn5/UjauO7qT7qCRuKMlEKNGSDITGa6YszuWBd5ewLLeAw9o14qaTu3NK75ZKNgkuZ1shj328nFdnrsUwfjKsEzec0JXGddNjHZrIdyjJRCgRk0yZklLn9a9yePTDbNbkFdKnTUNuOrkbw/u0VndKgsnZVsj/fbSc12atBeDirExuOLEb7RprpWSJT0oyEUrkJFOmuKSUt+as59GPslm5ZRc9WzXgxpO6cWa/1prWGueycwt4auoKXpuVQ4oZlxyVyfUndKWtkovEOSWZCCVDkilTUupMmreeh6csY/nmXbRrXIerh3XikqMyaVBbqwfEC3dn+so8xn6ygilLcklPS+GSrExuOLGr9niRhKEkE6FkSjJlSkqdKYs38dSnK5mxMo/6GWmMOCqTnwzrRPsm2sc9VvYVl/Luwo08NXUF83J20LReOlcO6ciVR3fUvS6ScJRkIpSMSSbcvJztPP3pSibN2wDASb1actmgDhzXowWpGrepFmvzCnn5yzW88mUOWwr20qV5PUYd25kLj2yv2WKSsOIuyZhZT+CVsKIuwO+BF4LyTsAq4GJ33xacMwYYBZQAN7n7e0H5QOA5oA7wDnCzu7uZZQTXGwhsBS5x91UVxZXsSabM+u27GTdtNa/OXMuWgn20a1yHS47K5OKsTFo30ioCVa24pJSPlm7mpemr+fjrzRhwUq9WXD64A8f3aKGJGZLw4i7JfOfLzVKBdcBg4EYgz93vN7M7gCbu/msz6wOMBwYBbYH/AD3cvcTMZgA3A9MIJZmH3X2ymd0A9Hf3n5nZCOB8d7+kolhqSpIps6+4lCmLN/HSjDVMXbaFFIPjerTg3AFtGd6nNfUytJLAoXJ35q/bwRtfrWPSvPVsKdhHywYZjBjUgRFHZWowX5JKvCeZ4cCd7j7MzJYCJ7j7BjNrA3zs7j2DVgzufl9wznvAHwi1dj5y915B+aXB+T8t+4y7f2FmacBGoIVXUNmalmTCrdlayCsz1/Dm7PWs276bOrVSObVPK84d0JbjerTQagIRWrllF5Pmrudfc9axYvMu0tNSOKV3S84b0I4Te7XUn6MkpXhfu2wEoVYKQCt33wAQJJqWQXk7Qi2VMjlBWVFwXL687Jy1wbWKzWwH0AzYEv7lZjYaGA3QoUOHKqpS4unQrC6/PK0Xt53ak1lrtvHWnHW8PW8DE+eup0HtNE7s2ZJT+7TihJ4tNDstTFmL5b2FG3l/4SaW5RYAMLhzU0Yf24UzDmuj3SilxotZkjGzdOBHwJiDfXQ/ZV5BeUXnfLfA/UngSQi1ZA4SR9JLSTGO6tSUozo15fdn9+XT7M28u2AjUxbnMnHuemqlGkO6NOOU3q0Y1q0ZXVvUr3ErC+zYXcQXy7cyddlmPlySy4Yde0hNMQZ1asplgzswvG9r3TgpEiaWLZkzgK/cfVPwepOZtQnrLssNynOAzLDz2gPrg/L2+ykPPycn6C5rBORFpxrJKT0thZN6teKkXq0oKXVmr9nGB4s28cGiTdw5cSEArRpmMKxrc4Z2a87Qrs2Scqxh974S5qzdzhfLtzA1ewtz126n1KFeeipDuzXn1lN7cErvVjSpp+VeRPYnlknmUr7tKgOYCIwE7g+e3worf8nM/kpo4L87MCMY+M83syHAdOAq4JFy1/oCuAj4sKLxGKlYaoqR1akpWZ2aMubM3qzZWshny7fwWfYWPv56M2/MXgdA64a1OaJD4+DRhH5tG1EnPXGm5paWOjnbdjN77TZmr9nOrNXbWLRhJyWlTorB4ZmN+fmJ3TimewsGZDYmPU1jLCIHE5OBfzOrS2jMpIu77wjKmgETgA7AGuDH7p4XvPdb4BqgGLjF3ScH5Vl8O4V5MvD/ginMtYFxwBGEWjAj3H1FRTHV5IH/yigtdZZszGfaiq3MWbud2Wu3sTZvNxBKTp2a1aVn6wb0aNWAXq0b0L1VAzKb1I3pL2h3Jzd/L6u3FvL1pnwWb9jJko35LN2YT8HeYgDqpqcyILMxR3ZowsCOTTiyYxONr4jsR1zPLosnSjJVZ0vBXuas2c68nO0s2ZjP15vyWZ1XSNmPmlmo1ZPZpC7tm9ahfZO6tGiQQfN66TSrn0Gz+uk0q5dO3fS0H5SMSkqdXfuK2bW3mK0F+9hcsJct+XvZXLCX3J17ydm2mzV5u1iTV8ieotJvzmtYO41ebRrSu3UDerVpyGHtGtGrdQOt9yYSASWZCCnJRNfufSUsy83n600FrM0rZO22QnLydrN2WyEbd+7hQD+GaSlGnVqp1E5PpU6tVMzAHRwPPTvsKSph177i7ySO8uqlp9KuSR06NK1Hx2Z16disLh2a1qV7qwa0bVS7xk1gEKkq8T6FWWqIOump9G/fmP7tv79VcFFJKdt27WNLwT7ydu1j6669bC3YR+G+YnYXlbB7Xym7i4rZva/km2mFZhaaQmhQp1Yq9TPSqJueRr2MVOqmp9G0XjotGmTQon4GzRuka5tqkRjRvzyJuVqpKbRsWJuWDbWsjUiyUaeziIhEjZKMiIhEjZKMiIhEjZKMiIhEjZKMiIhEjZKMiIhEjZKMiIhEjZKMiIhEjZaVCZjZZmB1JS7RnHKbotUANa3ONa2+oDrXFJWpc0d3b3GgN5VkqoiZzaxo/Z5kVNPqXNPqC6pzTRHNOqu7TEREokZJRkREokZJpuo8GesAYqCm1bmm1RdU55oianXWmIyIiESNWjIiIhI1SjIiIhI1SjKVZGanm9lSM8s2sztiHU9lmFmmmX1kZovNbKGZ3RyUNzWzD8xsWfDcJOycMUHdl5rZaWHlA81sfvDewxbH+xubWaqZzTazScHrZK9vYzN7zcyWBH/XR9eAOv8i+JleYGbjzax2stXZzJ4xs1wzWxBWVmV1NLMMM3slKJ9uZp0iCszd9TjEB5AKLAe6AOnAXKBPrOOqRH3aAEcGxw2Ar4E+wIPAHUH5HcADwXGfoM4ZQOfgzyI1eG8GcDSh3ZInA2fEun4V1PtW4CVgUvA62ev7PHBtcJwONE7mOgPtgJVAneD1BOAnyVZn4DjgSGBBWFmV1RG4AfhHcDwCeCWiuGL9B5PIj+Av4r2w12OAMbGOqwrr9xZwKrAUaBOUtQGW7q++wHvBn0kbYElY+aXAE7GuzwHq2B6YApwUlmSSub4Ng1+4Vq48mevcDlgLNCW05fwkYHgy1hnoVC7JVFkdyz4THKcRWiHADhaTussqp+yHt0xOUJbwgqbwEcB0oJW7bwAInlsGHztQ/dsFx+XL49FDwK+A0rCyZK5vF2Az8GzQRfiUmdUjievs7uuAPwNrgA3ADnd/nySuc5iqrOM357h7MbADaHawAJRkKmd//bEJPyfczOoDrwO3uPvOij66nzKvoDyumNnZQK67z4r0lP2UJUx9A2mEulQed/cjgF2EulEOJOHrHIxDnEuoW6gtUM/MrqjolP2UJVSdI3AodTyk+ivJVE4OkBn2uj2wPkaxVAkzq0Uowbzo7m8ExZvMrE3wfhsgNyg/UP1zguPy5fFmGPAjM1sFvAycZGb/JHnrC6FYc9x9evD6NUJJJ5nrfAqw0t03u3sR8AYwlOSuc5mqrOM355hZGtAIyDtYAEoylfMl0N3MOptZOqHBsIkxjumQBbNIngYWu/tfw96aCIwMjkcSGqspKx8RzDrpDHQHZgTN8nwzGxJc86qwc+KGu49x9/bu3onQ392H7n4FSVpfAHffCKw1s55B0cnAIpK4zoS6yYaYWd0g1pOBxSR3nctUZR3Dr3URoX8vB2/JxXqgKtEfwJmEZmEtB34b63gqWZdjCDV/5wFzgseZhPpdpwDLguemYef8Nqj7UsJm2gBZwILgvUeJYIAwxnU/gW8H/pO6vsAAYGbw9/wm0KQG1PkuYEkQ7zhCs6qSqs7AeEJjTkWEWh2jqrKOQG3gVSCb0Ay0LpHEpWVlREQkatRdJiIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIxLHzKy1mb1sZsvNbJGZvWNmPWIdl0iklGRE4lRwM9y/gI/dvau79wF+A7SKbWQikUuLdQAickAnAkXu/o+yAnefE8N4RH4wtWRE4lc/INLFO0XikpKMiIhEjZKMSPxaCAyMdRAilaEkIxK/PgQyzOy6sgIzO8rMjo9hTCI/iBbIFIljZtaW0O6dA4E9wCpCm8kti2VcIpFSkhERkahRd5mIiESNkoyIiESNkoyIiESNkoyIiESNkoyIiESNkoyIiESNkoyIiETN/weoNT2b9e586wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value of C is 3024\n"
     ]
    }
   ],
   "source": [
    "\n",
    "# Now we have the data we can easily estimate the hydraulic conductivity from \n",
    "# for values of C and use this to find the optimal value of C for each formation\n",
    "\n",
    "df_opt = df_coloc_bores.copy()\n",
    "\n",
    "\n",
    "# Do a linear transform\n",
    "\n",
    "df_opt['K'] = 10**df_opt['log_K']\n",
    "\n",
    "# Scan a number of different C values\n",
    "\n",
    "C_vals = np.arange(0,10000, 1)\n",
    "\n",
    "sum_of_squares = np.nan*np.ones(shape = C_vals.shape)\n",
    "\n",
    "for i, c in enumerate(C_vals):\n",
    "    sum_of_squares[i] = objective_function(df_opt, c, equation= 'SDR', loss_fn = 'L2')\n",
    "\n",
    "# Explore the results\n",
    "\n",
    "# Find the C that produces the minimum of the residuals\n",
    "\n",
    "best_C = C_vals[np.argmin(sum_of_squares)]\n",
    "plt.close()\n",
    "plt.plot(C_vals,sum_of_squares)\n",
    "plt.ylabel('Sum of squared residuals')\n",
    "plt.xlabel('C')\n",
    "plt.show()\n",
    "\n",
    "print(\"Optimal value of C is\" , best_C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATkAAAD8CAYAAAAMs9NCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9bbQtV1km+ryz1t7nnHwAwQQ6hoSgglfQIYYDKuF2Y+tJgsMe6LATkwbECxLGbbQFoi3QOi5DoEUIJ30dfrTxQmO3CB0ueOV6IeSIKI3aYEKjECISETUmnRBRCUnOXqtqvvfHnO+sd86aVavW3mvvXWef+WRUVq1ZX7Pq7HrW837NScyMgoKCgoMKs98dKCgoKNhNFJIrKCg40CgkV1BQcKBRSK6goOBAo5BcQUHBgUYhuYKCggONXSM5IrqCiD5LRHcS0at26zoFBQUFQ6DdyJMjogrAnwM4BuAuAH8M4Bpm/szaL1ZQUFAwgN1Scs8AcCczf56Z5wDeBeC5u3StgoKCgl7Mdum8FwD4G/X9LgDf2rfzueeeyxdffPEudaVgL3DPw3fiS186Gxv3PognXfI1+92dbeHPb/s8aFYBZAAigABAPuHacshZQ9xZGdhHfckd1znFGOsr01fKfKH2/r68de/9zHzeiJNncfl3nMl/96Vm1L63/enWB5n5iu1eaxXsFsnl/hqifxkiuhbAtQBw0UUX4dZbb92lrhTsFb7uptfhgv+0id//wE/ud1e2hWPmSlRnPxI4dAg0mwGzCjCO8HhWOUIw6k9bSE9IzmpiSogo851kf+Z2O7M7j26zNmyL3Et2CdmpvhIl/fb3hcoAxoArgw9+7s1/NXzCYfzdlxp8/IMXjdq3Ov9z5+7kWqtgt8zVuwBcqL4/DsDdegdmvpGZjzLz0fPO2/aPR8GEcM033Ib7nra5393YNk7Yd4MbC9Q1uK6B+QJY1EDTgOoGaBqg9ktj28XycsIhahf/nY1bItIRIvXk44jIk21VgYwBEbmlMm7fviVcOkPMQESstKz/I8AA7Mj/9hK7peT+GMATiegJAP4WwNUA/tUuXatgInjdN/0/WNjvxxXn/Rlu/uKv7Hd3tgU+uQVYC9powFUFmMqRyWzmyMK0hMRCJGZAK+RM3EQBskGs5ABHOiZucyrPxPulam8MpL9EjpzNisf3gMFY8DhzdS+xKyTHzDUR/QiADwKoALyNmW/fjWsVTAv/+qv+EC94+iv2uxvbBjcNqDFgtD4XZusJzsAzAkAEgnFEJ0STElr63aSER60CJH8e/8kVxYTnOgKQmLH+WDmuqvK+wT3GXqu0MdgtJQdmfj+A9+/W+QumiYsedw8eeNzx/e7G9mEbcFOBALAxPvZgwJVVjmbjuI4ZZL0SGyK4nB9PUGlV59fFJ9chPAI1VhEbxYTXhxz56X6ImtshGIxmAkSbolQ8FKwVX7n78ajm+92L7cP55Rrnm5NPdn43ttaZdbngAJAnEyG4xB/XWYyJgwGyGBNt48rEPrzw6Y/LnTu37BIseNSyl9g1JVdweuJffu578VW3fWm/u7Ez2Ma9hoacogNAVQOgUuYiugSnScwfH333nzxEMkTOz9an1pjBymTtBAyUmevuhbukliPjNQUemj0msDEoJFewFlz9R9fi43/09Xjsxy7EH33qx/e7OzvCCftuHKuuctFUUT6WAbJgJpC1zgdmHQmR5TYIkUIpp0BuRm3LIHCQddtDcIHYRXOFAJUpS5Z7lKRa1/6/XcJeq7QxKCRXEPD4t/8c/uqHxue4XfzrP4uzbzuMx9z2EGb/Y47PP3DdLvZub3GiuQnHZlcDjQUZ40xWFrLrUXMaJqPeNLktUXNOxbnzMwiw8MEOdJQd4P2CouyIuv1ibvu0BtWWAwNYTNAnV0iuIGAVggOALzz/1cDz2+/P+Zrr8HeXfjW2HmVwxn0NHvHZf8QHP/m6NfdyD2GV8VVVbbQVPkQgas74AESVEBdRltw6bX0I6SHwljKD2ROYBcgosm2cr7A3HUW3aZN2jYTH4EmaqyXwULAWfPO/uQF3/My5qF5wH75y6UO49xkG9z7zHHzXs16/313bNk7YdzsFB7iAgyT9snWVBxKEGIOU4FTQQJKC0yUEHCoCKnIkKsf575wGKzJBhpBwnAY/1g0GmpHLXqIouYId4ZKXHscZ//J/4qf+zYdwyaG7sQDhF774bHz0rK/BP87OAfEZuOzpr8Utf/za/e7q9sDsTVVXUpWlBh0o6JRPIU9wRGDZVfnzxG9HzK0mss5CBhEYtjVdfd4cM0CSgysqLVFobNCasrrPBmtTc67iYXooJFewbTz1t38K/+XfvR3fdNFdUfvRh78d933V2bj1K0fw0D+egQe+9mxc8aR/i5v//E3bvtYxc6WrPjAEUGuAcL3AieambZ93DNhyG2zo882lSJRSluAMtcGINAAayAiReUkw4bpsvY8ukB+1QQrvo9MExlLoIEy0dv8Zocn/DOwrirlasC1884/egI897Tc6BCeYkcXGRoP6DMbJRxEWFzwKz7no5Stf55i5EpcdeT7M4cMwRw6DjhwBHT4UFnPoEC479DwXJNgFjDZZh8guZxoKwXmy06ZlXNPqFvb7OxM1MVlDrp0jR9Y5dzlTdahfO4ALPNCoZS9RlFzByrj8qT+Nf/jRs/EVu4VDqt3+zyfiLxZfwe0PPwv3Pnw2FvMZyAJ2RlicNUN13qNwxfkvA+q6TbT1JMFNE8gD5ArPiQjmzDNBG7PQ1ikwJ3IF9I3BZZvXgOvaEdOawcze8c8gtujVB6nJKkhUXCA4bbZ21F98zjY9TgUfmJzJChPM10jlybWVHRlM17Cd1mKyujy56Sm5QnIFK4O2amzeN8Ob7/82POnBS3G2OYkFV/iH5gLcs3gUPn7/43HP3z8CzZc3cPikS2toNg2aMzdgHnk2aFEDdeNMQE92xH5ED1FNYpJWfvSN1GGuhh0iQ1EJ1rHqqvWasOLjIm+yRm3tJ1nuRlihTE9lpkYEp6Oten9IjpxLIWESU9OllZDxpjQIMCmxue8gn9eXJTrk0012ALvHKm0MCskVrIyb7/hZPOfrfgK//aVnYescRn0GgzfY13MC1VcMZicJZz0IbDzA2PyKU0HN4Qrm0WeCFg2otqDaKzmv6KjJmH4ZcpNPN1KHAVsLahpHiI0FLZyqu2X+zrXc7wn7bqcSmUBN06aQELm2Zeaq7Kv8cK052qq6nH+OhTR8CokjKk96zF64Mah2vwssfVEpJqnKQ+PYLhDdmlCUXMGBwgfufDP+6b94Ex46d4bFWQbNYcDO3ItTzQGzAKqTjNlJoNpiF/0zvvYS8P4n49Wcc+xz3Qz7t4DuiBzWvcBMXvFZdkMkrdmpzpZB1KPmwk4cL2mFA2IztUNwwuOaJ8hHVgGXFOx/SMgwYEk0nXuuQm4WsbLT67pfOiq8jmcEQjNBN38huYJt44zP/wNmD56N+SNnqI8YNJsEOwPAgGkAUzPMglEtGNRIZj5B3noidi44+QS6SapJfWhbxuTbZGw0ZscOnuCIaP1mq8D60iohPW269iF1+qNLcLmUEsApJHd+uJ0YISOZxI3JyoSF99XBOnIjZ86iQWsaNz2+wx2imKsFBwt/ey82G4vq5Jmoz9xAc9ig2TTgCk5tWEd01HhS8kP6MCO8yATj0h/IqyP/wrbypZvzBXhTS2WVEpFTc8YRHFcVqKpwzFy5nkCETx0R8zQKQFgbyJasL6DPQZOXcHOG4LI8UfnAg2U/aAAh6DjDvqY28dUJufkqCaZEzYWaWHfGnYJBmHO14/OsG4XkCraNm//hrbjisf8aFQBaHEZ1eAPNoQrNodZkIYZTVuISInLvFKm8rZD75T6Z0eZ8kVMkgP9I/XXGM4SsW/g2gzAhzTrAHJuskoKRmqiCZcrO3XgMowIRuS54i5aVqgvEJrdsE0UnAYtUzWH9Ss4lA0/PXJ1ejwpOKfADD4AeeBDmgZOoHpxj9tACs4cbVFsW1dzCLGxQc9EEVMH/1L7YId9L53zpfC6d89UHojYy6xOHj5kr13jDtnXuA12zeolPEeq+5XsbkIBKMRlYKgr7tyVg/lijTGCdcqOrLuTZSPsa0fiE4GXLXqKQXMGO8MGH/gv4wYdADz4MevAkzFfmqB5aoDpZw2w1MHMLU1tvsnKbkQ9E5MXeVZcSnc4vi9Apn3ITvIQcu8pP+OLnadgp0Z2w73YF+5IMbK3L7VMDbEqRfCiIV5HOdqDN9pyc3n8gfLT3niySDMzGJwJXBDYI3xGShdvkYSQJw9nk4DWAmdCwGbXsJQrJFewY8rJTY0HWOjJTC6zzn5F64cWMjRBFHNfwEhpPeGkS8TYRqh+UmmOpbQW6pmvfGG/LEKnceIGY+poUVfWEPDsO+ylyA2I1p7GmCKsFjVr2EoXkCnaOKMLplmCeamWT+tR2E94fF9TdGs1WreZcA7clXpmoMI2856DiEJuzUO26LVV/ywbl7B2ReI0pJHOejVr2EoXkCnaEY9VVcYPVak0RnVJuWRW3DXCq0FI/VIp1mGdpv3U9q2xPid23L4tDhEtkkoLle0RUQnRQ0dmsmsOupIukkMDDmGUvUUiuYEegys9L2ilFAqB9cCrC2jrp/b52PAH0dyR5ifV4bMFk3bmay5msACDDMUWBB/0p1WpCgsyjFV6EQGK6jeJtwPg3O+3rDtEwjVr2EiWFpGBnqKp2hndNNMwgJudnD++P9rqvYsZlis31sYZccqtR1Q/c+Nw5Cjl0ZGxbJrUTpOkkIYfOjzenJoZ2OXXtcV151gP9KBVhSRK0eybuxAy4+l35sSBpTMjQwMu6VI2uz1wtFQ8FBwrHqqtgzjgDMI7kJKoXwbLL2SKK3i2t9gJWUROSo5Z+Am3uHODIV9RVA+ix6NYGy24+BvHTpflzoFCR0Jk9i+V4n+xLqp1igmvvHR2eGtxmKLgN+u9hPUWsdo8jp2MwvR4VnDIwZ5wBOnIYdGgTmPWpOR9ZbWy0iB+rY8JqqAhrnFqi1tM2GbHEK0ySUUyqCjSbuSqI6qodma2pycrMLoWErZub1RfBU2P9p5Cs3K9ft2K2qmeQKNyxBfSpQF0aZAgRYDc4Aq+B5BhAAzNq2UsUJVewLRwzV6I65xzQxgZQmXbSY3m5fJE4AISx0NJUt7HCLRkmqD0B9X4SfGDCGlcb63PpnOrrO+EK0CartW4EFEkbkTKvhFCkHKxz6VTN6XbyFvEYXhCzvae/3euKmrZrMVkZhMUEy7qKkitYGceqq2AOHwY2N4DZrPXJVWYwB0sUS1Ausi2j4jglLmlLlFsUYQ3DF5l4P9O2EakgRBoZ3gkyEdZo3UKRnb9/i8Bo8TPp+ivJjlB1fUSVtEuwJ2zTuX47gBOoJRm44BTHsdnVMEeOuGHINzaAWeVUnCi59C8qfdn7aj2HkDuv2hbNRhURnZCdcWarMa0Zuy6iU5PchO/p/eoos423RVHWQDhoiS73iHYiuvQzjxKa1xF8GJcIXJKBCyaLyw49D+ZM5YfbmIGVLy6qM9UY8xINqTlBp+wr45sDwrpLbZFh033Zl5R/JUS3lkRhT3BRYnDGZG19b2jVXKpwPUiRXlikPdovc1zUt7jiJOTyybDzTYOdgrE+JUdEFxLRh4noDiK6nYh+zLe/loj+log+6ZfvXnau4pMrWIrLNq8BHToEc9aZwGwG2vQKbqbMVGPi4YMA75dTJ9qJWlDR0yilxLezIT9YpPe3ydDehuBCnw2Ayo3s6wfr5KoCoQFbd8zKwzKxn/U5auMwpJQQCcsk0BaOjawM14SW7Bgh+kqGWrFGPg0ll4eojw/XT/qik5D7VLX19bdrwBqDCjWA65j5E0R0NoDbiOiE33YDM18/9kSF5AqyuOzQ89xsWGeegeqfPDbMtcDid6uMG4E3KaDvqK8x7076Xqi4QERogOMU1qPeEmSQOvbDexPg3nyuQMa4InpDTq0YH5RgdkOmC3n6SXSOVVcBzINkd8xcGcxgd49q6KWQO6dSV/x2IuuGIocFmvamfXyh/Z8MpSRVC4ALcKQIBCkmsSPAUC8c1B+HCC9Z20Z+6wZc18CiBs/ny/6VloJBaxs0k5nvAXCPX3+AiO4AcMF2zlVIriCCM0mPwDzqkU6xHdp0PjdfOQDjSW4Zua2CMMqvakuJrpMLBzXqLWIyTYMWlVdwVeWITn8nP2+DAWBNsCGF7ABEhDfWrGUhlETNtQqOg9okBrhhN24ow4vOdmhzIO+GS8fq0+Zs1zyOTVWoVBf2s6ftFAxgMb4u9VwiulV9v5GZb8ztSEQXA/gWAB8DcCmAHyGiHwRwK5za+/uhCxWSKwAAXH72D4GOHEb1mHODr42rCrwx6ye0IctEkU5fZUOv6hs6rxzDYhKqvggpiunqJ7YBEBMbM7TpGvrYlhO4Inx/wg6xSTAjmeg6f5PKRAyzexGkIsSpUviKBSAk83qz1SlUgHLO+sgUVc9ZqzhdS6xVHDs/HNc12H/uHCuNFXc/Mx9dekaiswC8B8DLmfnLRPTLAF4Hd3evA/AWAC8aOkchuQIAAJ3/GEduxoBnysdmRvhYMk7uUddUPrYscmouHCwEkCn5Mm5IdWLKE51RPjq2fg4FldtmjVdVpEbh1aydPJNlzygXUfb+OXLdaP1zDdz1rCI441Rd9/khpIJIgEIrt5BobG1XxTU+AbhpwAs3Vy3tRI1DOHV9sUwi2oAjuHcw83sBgJnvVdt/FcBvLztPIbkCAEB93iOCHwhAf+KuSmsI+Vbe5Arbo+MyBJi8TKljPQpYDJmtJLKn7XfYX0dbE6ITJQc0LaGRqj0lX3wfyM764IS/VBhV13TL2OQec3OyGvVMWNW0arMVDLCavwEI8zpkS+ZSJddwrOJCJNVGVRhoGqCuHcHV7+rewzaxrlF/yTHuWwHcwczHVfv53l8HAN8H4NPLzrUjkiOiLwB4AO5nsWbmo0T0aAD/FcDFAL4A4KplNnPB/mP+qM2IqGRdxoWLlEIY6Rb+5WRv3nmzSkzJvgx7TogJS4hOnyIlOvHNafI07n/OL+YPVERHXum5yUz9JAho/KpPMwGyZOdP4K+j+jtW8Qr5sVNxXKFjtoLDRIPKK4fwjN3zkka0uXbBTG1VHBofdPBz24YyrrpxAYfFzgMO7e3ROpXcpQBeAOBTRPRJ3/YaANcQ0VPhaPwLAF667ETrUHLfwcz3q++vAvAhZn4jEb3Kf//JNVynYBcxf2QVZtVyi1+vxeTxETzrZoci66OCauKUkEWvAwSCnAmbIbsIQnQJYUZEByDMK5oSHXuzVROdKBsAIAY3/kJV5RQX0Oa5pWSXposAflDOAfWSmqods9VHZC385NAAGoAqF6109yT3DLSM568vyjrJf6PG3QPJvdgkolrX4Lpez0xmqkvrKuti5o8CWVn4/lXPtRvm6nMBPNuv/xqA30MhucnjoccYmAXDLJxfyNSAady8qaaBL7JHO1dDwzDkzCOt6tDY4ZrU1K+m2sL7m1N0WiDYdp+wp1ihSl3KNUKeGruBAjCDi6Qyg6hxZCgqy5Af8LKt59SR0g50gvHQ5DASZfUmLDXsVJy17Qz3bNw1Z8b/gnjSNnCqLmPmy/NwwQT4Iei92do0juwWNWhRA/MFeLEAb83B8zluefjXc/9COwDtecnWGOyU5BjALeSSeH7Fh4AfKzYzM99DRI/JHUhE1wK4FgAuuuiiHXajYKeozwDMnGBmgFm4F5BreWn9uyZSQrhJrDdxlrMjiU4NZWqmpu3LfHQphvLo/KAAfeklXHl1JgnDMhm1dbYeA60J62dudmlvKicuup4iuD4fnS++D/ccAgAEGJ9W4hUdgfxjaQ1WWE/oJnl+EjX1PriI4Kwf/aRuQHXjFVzjiG4+B29t9T/fbcK5Adfjk1sndkpylzLz3Z7IThDRn4090BPijQBw9OjRod/+gj1Acxh+ijv42aAANs48ZWIYcuRmwEBNSmlAcnGjOEAHOqooMGpbFEwYATnWZsxXAGFSZUnQ9TlhAIL6itqN8SYqOZUH66/hzVhRdlXlu+zVp1ZwUf96AhLis/R+M/bf5TcCQNsuHjkfZW3nTG13Dj64lOBq739rvO9tvgAWTsHxfL5WM1XjwA2aycx3+8/7iOg3ATwDwL0SASGi8wHct4Z+Fuwy6sOMihy5tVPgEZrGE1wYxdERnSsooNBM5De3Qm85dF5cQnRL1ZzA81HHfIVz6sNqk9aZgwxv3iXtMkQTyM1IL1UQkERhFXzo5K2FGllCNndO7suye27S5tUciBCGp5KBRoUQqb2epJKk0W4dcID3x6F2M6i5qoY5+OQW7NYWTjQ3LX+u28A6Kx7WiW3TLhGd6WvKQERnArgMLpz7PgAv9Lu9EMBv7bSTBbsPe5jRHGLYQ4xmE2GxG4Cdwc/xiXjeT09sg6ymTVQRUvIyJu2pWTs4PLp2+Ctijiak9uZqdoimdAy8ZJim2ASt2iGbyOSXVaY9FP+f3GMwO7mNkIYBNX1gIag/JPO6SvBBBRskXcSbqc5EXewqwQmmOJHNTpTcYwH8pk8gnAH4DWa+mYj+GMBNRPRiAH8NYI3TlxfsFuxhX2/ZThAAAOCZyx+zDIAp5Hux5JaR5HChkzKiJ7GJJnJJtnNQMcibrmk6iSYTvS6pK0Bb2wq46CsoFMZrnyFXpiUNMV+lCkL8dBJZZet8a30BCCAQXjbqmprU2iepgi8uX4+ARpKBvarTx+hn6M3U8Jx1ukjt6lLtyd0nOGZgYfeWwMZg2yTHzJ8H8M2Z9r8D8J076VTBPuBQ44nMeD8QAHYqDtaTm39vyHqfHXEoR3LRv2FRF7Lxo8Y2EsqWYqJDJgiREodOV0l8coHoQk6db5NUEm+yskGcZmIV2VnrRxkGgAqheD43yGRf4EFD++Us/GjAPvAggQX9QwA/QQ04mwwMIAyjHgIoEkVuXD4cb22tNeG399aw1jy5taFUPBQAAA6fNcdWtQFrKvDMwFTGDfhbC3OJw1xSLHzamfUcQu3sUQCSMiJFcAMRVgKAxis7HUnMmYG6MqMTnUUIOITLhTw6blM2hPx0AX3jrx36bpQpaVsS9snG3W6JoqPeBOGQSgI4ovURXyE68goOgPsegg2Z5yZkWXsztW5cusiiBm/Nga0t3DJ/Z7Yfu4F1VTysE4XkCgAAmxs1mtpg0RDYErhhcCNRVjjflrdm5cfaTfcHdEKq2ueWElzfEN3CB5Qomz6TdQBheCL9wmllp9fl2pKSIsX9mW1I0kh6X+dVa0C1qRpSSfymVOFp5aqeLYkyDaOLNIBdV+H9yNvAwUwhKTggOLK5wKKu0NQGtjLgGYMXPtigCU7+hjP5Ip0kYOYuwfUFEzSJLcubS+tre3LTOqN3CLl5pZglOh9h7iU6oEN2Wcg+fYSXBl5k92RTpG7dga5dKzmJpuqqBim+X8MQSuNRzNWCCWOzajCbNVhUM9gZg2sGV9zmzAnB9URT5aXr1LamBDcUMU1f9lTN5dDnA/NJyZLTq0vA2uGZMkTnt0VER4SQPBxdY8kLPUhwmW1BoSXmerp3YqqSShcJuXG++J4Xe6fkAOz5/A1jUEiuAABwaFZjo2qwVVlQZd2QSxJsHSA3SWMI69oJniM4IUPrCLSDlOhGvjS5nDpJV2tnlm8VYBR51USnztMhOpWA3O3AgA9xKBihTE/ZL5jrvn+d8/v9w/MNo4uoyoaFr03dg4CDwEVXpzclYSG5AgDAhmmwUVlUlUVtnGM8MlEH1BuA+MWP1FxMcKRMVmpalRZbOepiFgg1nD3oSxoOlRCGWqID2pFLZMfIX0dxGVgYHkmbvck9p6khwFK/XBh6SVcwaIKXczbx/UVmqvbDCcEllQ17iakmAxeSKwAAHK4WmBlHciTJtSQmKEJKSbt4U0kK9mVpbCgt0o7wQG6pHyukirgEXQDuOKI26oieYAQUAWRVZrudDPwAmDLhjdSLqqgrAwRV7qUIOuTTASoi3KOyUiwbpUSjSQjTz0nRmdNV0kRkAMyTJ52JenILdr7YUwWnUczVgkniFf/jB2Do0TDE6EyYwu1nOjl0u8hLx91oqhCckFsaXVUKTg8sGSUEp0opp5xySGMjMniA9s9pziVKTFnkVd1AvGG4P0lARUz2NGqq0ag2TW5KxXFdh8L7vahq6EOJrhZMGsaTG8sfKcNHG9xLrdVcmM1dCMwPrJnzw0UlSLnoqviigr9OmW+Sm4cM0Y2FJzptugbiSKOszBlTFh0fnPPVDfRnyD+Xdk8TXR+UH1OUG5idcmtcVYN9+CR4sXuF92NRoqsFk0QDE/0Cs5AZJ+SW+R7a5cXrW1JTD2jVi3a6a1+VHqJoLLJlU4iILqSWyHZNbL4tO3eEbJeBDOSSqyq7PkWX9LtjHjO7WtSmTRFxZVsL3LL1jhU7sX4wE+pCcgVThGXnMG6YYC21Ks4iUm8uXUHN79noSCo6wYZIxQFdtZIGLnSAQZtwFn7EIx4dbQUQm4Cp6SrbVT+Y2hy6MIST7JtRdOE2kvd6O+YspWZ8+szEPFXqzY0usoCdL1a/4C6hmKsFk0TNFeZ2hsYaWGvAjYnMUvJkF5Sb5UB+LRG2o2hkVVxYTxggzTXryyHbLhJfFxOFIEQnf852v48lOg3uK+KXLvVVfQB5cpPnJ1MH+uGTwthw++SDS1F8cgWTxcPNBhZNhUVj0PiyLmoIZJVPTuZ8EKJrVMBBBxkyqSIRwQVllTErPcKYaznC86puJYTgBSI1FxKFZR/D8cCUmuhEWQZOpv4SNX1NDVGMyjTPbe/7cWAZVcTnwE2J4ASF5Aomib/fOgNfnh/Cw1ubqOczYG5Ac4JZwC01/DwPfrGtmgtzf+YUnH5RUwUnpCHzneq2HPS2vv2GnPyK6MLEO/CKThFpNNCmL/IHu6J+QOWpuS/5fmpov1pqHodorWZe9dxkZi3xvzV231NEhlDy5Aomi4frDWwtZqgXFXhhQDXB1BQIDeKP49YfByD2x00MoUhf902pq2yicBiSCZnZv/xxtn2Js69zquD6lJ/yObJMniNIVbA2Vxs7SYITlDy5gkniwesibkIAACAASURBVMUm5vUMTWOAmkBCcI0QHQe/nFMa6DrxhyDDDg0VtO8EQ7m2KdnJwAK8DaIDoppSTXhxfyhWfOFY/5lWSvT576yfCtHPL8HMe1xwvxqYgfogDZpZcHCwVc+wWFSwCwM0BKq9WVqr4APHn6OQGakk3oYw3HjUtkPkLKbI/wZ08+dkH52nZ1oiC+klPYTXuR4oBDI6hNeTkjL4XGXE0txgnRNCMVcLJomtukJTV0BtQNpcVdHVKCEYeXeU20DtkibfpmpOR1ajwvglpLdNMuwjOtcXZU5STHSdKOtYJCZudHxKdGOxyr57jOKTK5gsFosZbE3OVJXAQhJgGCI3mYfAzQ8qyoRcUXtH/aRJZRQrue0SWzrGXA86icDoMVvlnAnRhcsNBUl0n5j7yS5TNkaSGN1Ml8yGwIXkCqaI+ckZ+OEZzJZBNSdUc2eymkVMdtqccuPLkSsgNwRYchMwVwzyfpmQACIF70CmdjUmONazZoV1uaZe98fqTBT9fqUElCqgZBBKSTIOOXRQ/WL20xEqYhyr64TYhCzZD8muhmPXwQWGdak7/hm4Cg2ZEcwsJ9Z9Rgk8FEwSdlG5IIP44oTYuF2yEKsuqhRQ5iq5kUWiOU5Tv7RScYHgUmW3TO0t+y5t2ve2RChFAYvcIAFj0aviMgN2yggk/scjTI4NoJ0BbHqOfYHj6umR3HSfWMHeoSZARVMhVQ6qrCuHdj5TuORdTUiKsLhSKmTsYqi9BmIVB8CZltqPJ6vLcuX62vQm03PenSxajao2XnpcPAMYLZsNbF9BaKwZtSw9E9GFRPRhIrqDiG4noh/z7Y8mohNE9Dn/ec6ycxWSK/ARVQk0UFzGlUEYMRjw4861EzpzpSZ3VpM3RypNv8CeBKPtOUUnn0Qds3QUwSHep09wBOLMEV12/+El6rtMCJS0cXLfrcnuPinabnDMTHcqY2YatYxADeA6Zv4GAN8G4GVE9GQArwLwIWZ+IoAP+e+DKOZqgc+Lo7YIXwguUXHy4hIAGQ7dvaDsLDqCeyErctP+WeuITuY3HZqYOacANQkkvjgm6Yg6hyD96V6WdRHM0fh+A9FZHiS6aN9wzUxaCtAGMjjzLOQ+E5MfYqbKvK4T9csx1meuMvM9AO7x6w8Q0R0ALgDwXADP9rv9GoDfA/CTQ+cqJFfQJv1KhYPKicsf4IgnuIyEcAz5aAPFE8boCZs751KqRkYGVgTnzDW0L3uCzryrOdtkZGF9e84kt02nl8j3ZUgIMhoMgDkp/OfuUE8SnQb8M7DOL1dVQDW9eRQAuLzDXQgKE9HFAL4FwMcAPNYTIJj5HiJ6zLLjC8kV+GL71lSlAT+cEJzM6M4WIOMjhl7ViYnFlTs3AE8+1I2uAlGENZCbVnPaTNUqThPfMsdLbvBLIR9/7TStBEBLdgPE1jfHRBikM0d0Q0nSui8+AEHRs5iul2mF6Oq5RHSr+n4jM9+Y7kREZwF4D4CXM/OXqedZD6GQXAHMwhXjhwqHBgjDKGkEk9SvG6fYmEWFwJup7Sz1HFIk0KZRJAgkUbWkpcnOXScmt6DgTHKODFqiQjucUi7SSom9ysPn7V5I35QiytysYZ3oKsdlZNpvadmpN+YQiDhWXTW5EUjYBx5G4n5mPjq0AxFtwBHcO5j5vb75XiI636u48wHct+xC0/1JKNgzODO1rXDQQYfB4IOsa2e+Ip9cJLGzVJ6oqla5ceX8eiEiO0Bwch7X2WRJ+geg/YtXbb1upOQ8g/uk+/Uct4rLKlR+6OiqBCImitxANLllGcjd5FsB3MHMx9Wm9wF4oV9/IYDfWnauouQK/KxbCIGGdlhz99fYS3QEnwjsfUyi8MgHIFK1kiTUptFTrtR3bZ4qZ3uq4Nz+PfelhFnHzybXSRWd8pstPf8y+Otnr53bXY88AnTGmyOilcZF2A+sseLhUgAvAPApIvqkb3sNgDcCuImIXgzgrwEsDTUXkiuAaaiNqqoI6+jRRtK/awlAAF2iC8doJZUhN38eTtu0SlTXzk4unURM2/IttKb4oOmK3fGk55C6BtKhlgA3Colfpgin0tZDcsz8UfT/vHznKucqJFcQkn8jU3UoupoDwZGP8AqhTZ2QOU4h23rMx1xwAV31Fm/rf6naEX2xPaKTbcD2yC6Y/O2x2Weqz63tuWjwUdsGbfoi1RPAFCseCskVhOgqIqLjxHRVn0sQoq/hAhSXN0WmZqvictHT1P/m9u0quXxHEqLzbQEp0QE+vcM35cguOv8KJmg6YKaenxZqPZ0ESA+Bzn6EYLb7PvVgH6bIvYXkCrpDKqUT12jCg29bpvRkRA89Z6ooKE0YytfGIU8OeYJL1Rv1O/Kdfw1dnxi5KGAn4pqBPnf2XodUXqapV8XpIZeE4PT8qn5EYDQWLJPZTBAMgi2DZhZMEaZ2iwyUaRoOycFhXU9mE4ZfEqUBSMrJIBnkIpq6RtW0hJeSW8d8VcfmwABgWQk4by57VadHESE9+GVutN8+ddd74VjBtT8OHEiN/KQ+1Pg2ITQ1n0OYX7WuAT95zS0P//qSDuwvJijklqeQENHbiOg+Ivq0austkiWiVxPRnUT0WSK6fLc6XrA+hKL8RhNZG4hwn17NhVm7ZL4HWc//efeN0gsgmKYtkQ0TXNhuYkLM1oym18qowHYflYbi60vjB5Qh59CgzVC/+zKCU+up3428YssRHM/n3Yc5JfjAw5pqV9eGMdry7QCuSNqyRbK+gPZqAE/xx/wSEU20BqVA4BRbS2CdQISOvCrVlqq4yL/Ug6wKE2IDsuZpZJYm5Nat8VRqUCtDuQ4SokvILiBHdkPIEJxr516CEzWHxg1tHvxwOYLzn5MHj1z2EEv/GZn5IwC+lDQ/F644Fv7ze1X7u5h5i5n/EsCdAJ6xpr4W7BJSYmsXjj+ZW7WXUXEkJiwQXvKseadN1CUEJzWyclyH3NBVbysRnWrrtAPZ5OEOegguPI+U4JjzBOf9cFmCW8wnG2zQOFWVXA5RkSwAKZK9AMDfqP3u8m0dENG1RHQrEd36xS9+cZvdKFgHur63rv+NLKsJpeFfVsQqTtp6CI4HyMQRUtweqTnTHVIpZ5rm9ukQXa4PY4huAL0EZwFq/A+BtUDj/XBilmqCqxtHbqcqwQGwlkYte4l1h0Jyvc+KU2a+kZmPMvPR8847b83dKFgFkR8uIjok7WJmuRfVKTqvTDzhLauSyJmpnR92Gkdw7TnUgsSUTa7Laj/3fSTR5ZCa5RkF15qntqPmWh9c4whOzNRTkOAA+B+9Hidpn9N0D7Dd6GpfkexdAC5U+z0OwN076WDB7qOat0RmGk94/tPUEmnlltxE2dU2vNCOAP1LnhtphMhNZ9DXiSgA0QYagJiwUrNzKE+OAbiMPd8vIToVdY0K6ZN8us5IJCppWF+EfJQ5a54KuWnzlBlUN22aSF27FBHbAIsa9qGHwHV96pCbwhTz5Lar5PqKZN8H4GoiOkRETwDwRAAf31kXC3YbfQTXklyG4OSlti3BRfNBrOuvPTsUefIZ3UyubYAJR6i6LPT99REcy7NLCK4R07RpCa6uwfMFeD4/ZQkOwCQDD0uVHBG9E24kznOJ6C4A/wd6imSZ+XYiugnAZ+CGL34ZM093yu8CAMMEF5unPtBg2U26IuZq+pILUgW03dEzcuSTCSDohN+2jKsl3r71zojAY5E7JjVJkxy44H/z+XBc1y7YUNfAYnFqExz2PqgwBktJjpmv6dmULZJl5jcAeMNOOlWwt+iYqCEIwTHBCakpf1xIZgWCiolPjpjgLHeTeKUuNfXFKYx5d/rMznCNRF126lT9OXpTYDLtki+YqzkNAQbtfxPS6yG4W+bvXH6jU8YEzdVS8VDQS3CIlF1rtooPLhCcIjdNEEzk2oXogNFqTkdFux1W5x95rkjBpcOZyzm3+4IuS/RNCc7aA0twvMeR0zEoJFeQV3D+0wixSTqJMlHFqQ6ol1uf1392ivJXQJz/NvaYHjU2Us1F10u5MKi2NuAAoD/RdyTBnRKJvqNQSK5ggogILqfcgoIDqGkjqulIGhGBpH44Cz8SSc4W7WkfQDZIkCOrPoW2TjWXmqnKF9dWMWiCazoEd+r64RIUc7Vgiqi2uM2Vq21McrVtfW9prlfj/qJzqkmSe2U7G+PXXQoHqvGkNqjm0m3BLKXEdO5XbOm2cKzePyIylTYSVK5P9E0qGbCo2yCDj6KiacAnt8D1YnLzNOwYheQKpohAaBkVR0qNpD64nInqTuhIwpmpnHz3uzCP9qmtFdpkzam5oUNzpitUuoiOpkZmatPmwTWNGy6paQ4mwZ2K0dWCgw8hOJOYp7JAlyWlKRJASxRRNUOX6AYhAYreTu7kDmMM+uFGIPLFpYGGdNDLpgEzxwRX1+77AcQUk4ELyRX0E1zDkQkWKTmgOwx3w206CBATm/jntJ9udVfcMLRLcLfettSE9eTWyY3T+XFN4yb7kRKupjk4PrgUJbpaMEX0Epy1ra9Jv7w2ITuBkJhSZc5N1uav8QSjbyniORnE74jofuMRV9BVcf6HQUb0FXLjxjqf5AHFThTybqGQXEFLcE3GB6ed6Sm59aml1DS1cNMRyubEp78W7ObLlSY425gEOypOJp0RM5XdsOWO+OyBNVX3o2RrDArJFYAWNhTca/8bNezqK9ULTLaH4MhpNjaJmtPzIAz55aLJqtGOPJJJFYkCFrm4h+5bCA7I9/xbOMpPx6rCQXIEc8GGNF1EjSpyon7Xkoucytj7EUbGoJBcQVS9EBFcolBI++CGcuI0lpFbpzM7vJdd88UlQRa5Vu75pH64oIJTSXgAUZRcwRQR/G9NQnDygqYvcApFZGQZvO3qhuEyrtVOtr0+jDu3jqYiO9Q5syc1aw9uykgOE+TxQnIF3Zd2SKH05MWF86RENTRTF6GzbSfoqLjUVM1hhVy5wcJ9HZTRPrnTTcVN0Fyd3iSJA7j8rBfiu571+v3uxsGDcp5HCa1jCG5VhMEv5bOnS6k/buiU0u9w8E46mFGUmahq+rw6z4ltS3hs1/PsTgG0c38ML3uJU4rkzD95DBaP2MQlLzm+3105MDj2zNd1VVxupin9kqpRgCMQucCD+g4g/JW10/7l1J0ivaGBMtNL9qi3rIpbQjSDImRpUKJ9ZhzWLZgZbPng5sWl4JHLHmISJPfnt31+1H7NOWeiPmLAFfCUn7xhl3t18PFdl76+HVVEqzignxA6Re0ZZogmio4JLUtaJtO2HQwR3BAGJqmOkN67Nu1lu/VmaTBdTxNTdcKYhE+ONma44lEvhn34JG7Zekd2n8vPeiHmz/wGNJsEskC1BTzlVTeg2QT+7Gdescc9PnXxz7/zZ326CLfpECn6/G7MndIt7XNjQy3BEbmifAPAGEd2ldsm+zEBqOLJpXXqSJh0JnXzLVFvbp/hexKyDfuJby4diSTlJ614/feOqSqJv95sPS0CDh7rMkWJ6G0AvgfAfcz8jb7ttQBeAkCm93sNM79/2bkmoeRABnTGEZgzj2Q3HzNXwjzibPDMj2RhAaoBs3BLwTg8+7I3upVlwmIoeNCzBILTi8+T43QdUEQmBEddghuDbRBctGkH6jH2A2Z8lvZgVzd0wHBlXWOW5Xg7upPaA8ANzPxUvywlOGAqJGcIOHIYdOQILj/jBThWXRVvPvts4IwjLmHTMkwNVAuGWbiZpp7yqmK6LsM//Z43tV9SH5lGjtDS7wm5tQrOAJVx7VVOzQFhMmlV49o35Plo87WP4NYVLMleU5mo6hocRVRPQ1N1TT65nkntt4VJmKtsCHzkkPv7P/MMmKrCZYeeB7AFzWagRz4CfGgDxIBp3AQrbAhM7H6KT6Mfy+3gf33umzt8wZ5XRDmRV08EeNPN72eceyCHjnoTM7TKqTkoM5baiaQ98em5UqO5UZcFA4YILged7qKfBfv+NclxBoBMqh2ukzn36aTYBrCCuXouEd2qvt/IzDeOOO5HiOgHAdwK4Dpm/vtlB0yC5EAEe2jDvQ+LGjAVaGMObixocwN0+BBsVYXMfFM7gmNfSgQmfNMrbsCnbii+OY1nXnm9G/VXtQmxELz/yRMOmADjhkYiISeLlrRyZVxa6Wk/nFZwgdS86iNEBMcpAWpSW8GUXIng9Lrfd0xZV1qk34Hkx2mcJqkjAeNv935mPrri2X8ZwOv8VV4H4C0AXrTsoMmYq/bwDPbQBvjwJujwIdDhw6AzjoA2N8EbM6AiUM0wC7dUC6CaM6q5fDKe+rKSWiL49h+4Pm4IRCLfhVTiTxgEonIEKOsmXoicaWrk0wRii0xURXDwwYeI4AxaX5zuq/rOypfXNqrdVyW4gXZOrj2IgQg0s4q0nk7YxRQSZr6XmRtmtgB+FcAzxhw3CSXHBNhDbpgKM990w09XxhU7y4sGwDQWXJOLPPg57Nw8oDvwHq8JT3/h8d6MeO1z0i9R3wt12//1ym3341uf9xaXcEmtNGGv1JiceiNRU+zXiYAKYBh/D7ZdpwG/lv930YEFITJRZVyZ4IdzSs4TnGkV3FA0ddnowTsiOL19lxTXaZMfB+x6oi8Rnc/M9/iv3wfg02OOmwbJVYT6zApmRgAOw5ycwczrdnx8AFRbmIdroGGYDQOqGTwj2A2CnRFM7f6Qv/X5bwFZ4L//xnW72uejL3KqMSiQGeBIt/+YXI5YLrr3tB92516F7J7+Q8cByfBgRaBGLsjgxl3MWaEEIgYadt8tuxF8rT84Le/S0GkjmuDgyUsTGyFWbyPIbemw6GlEdYigckot3YV9P9Rod905ItyPxio/p3y6marA2gbN7JnU/tlE9FS4f/kvAHjpmHNNguQAoNk07uXcMuANA7YGVFP7K2staGFhDIXhySybIH+jMiBmPPOq6/GHN/342vv5tBd7kzhEBpMdOLOeWlo5vxN3tz/th4+PIrqj/9tx9XLC+9zglRW3BGKgX2MgkpR+f+NHJCH/yqeOeI2qPa4lLoqHSTKIFVyVITf/HLIzcCXPZrBtCVZOF6HMdUg1DpHx6WiqYn1KrmdS+7du51yTIDkmwG4QiAl204BsBV5YNwmwH62WAVDTAIvWJ+6NWHcO415oneR56fdfjz94z3qI7uiLjoeXsmNmKt7goZdz6CXLHMcEXPKS4/jEr+aJTgiXSFGW9FGta1VFzpYFV3IhAlXk1kS4iYqx7MxeIFZLOtfNKH+ZMks75Cb7av9beJ7q+HVCkdDQD0sUdMhFWDWxFQxjgo9pEiQHAuwMABOqDQNT26ASZChpYgYW7YiqEdFZ//JU3R/XnRDd03/oONjIC+ra2LkC3Xr6Uibfl/2qjR2oMSU6qd0N/AMlOrRZRe13NvDEpu1ZdzSzUsyNcsMZcn47Vqov9J3a88uzyCm31DRVJVwpuaVDNPWlrqyK3ohtyl3+GSz7d2GiYT6W3LjTLa1kl31y28UkSI4JsDMnP7gCbGVQEbW1gI3/o6mbTsmjm7PYgGsOKsT5ftw+BMal33+9fAkqTJuy3/av3tJRabbyL2zSrl9S6XsOxP3b0nvvRAozx11y7XF84sZX4pJrj7fXVl3pVXLeJ+fGeWv3DPMtEMHCzbUK75ty7gH36YYq73Yocg8kKSHLSrSiQMxAfJ/NzoluaUqK4vzoBQ3t5NJtCsZhgo9qEiQnSo6sI7sw77AM+VPXgDGguvHvYfvXGtK5av9GGw4JrNELlPyBf/sPXO9N3JjIpK1DbMrE6hsZQzurB83W7u2rc8T7622XvPR49kWNlRwAQ+3glUHuevKCf5sb9/Iyw9eUOgHiVByBjH/zR0WM9TraZGAgS27pD0UvWBHdqi/PUM5dK2K3j2LCZrEu9b1OTILkmIDmkHsJZid98qm1jtTmc0dyAMi6CgjUDbAxAzUMri1ow9mp1sJFXCuEF5ur+DqazIS8NLm1Sg5dspP9gIj0onsZenN7HOjOH8bZ/XKHxGkTaCeGEbazQrhuZ2aAjPtxIMsgJpBPKyFmUC0zxlO4YDtuWs+tKCWX9bGlSk3aZb8+hOtz15zsw1AAoOffKbSpX4goaAP1oyU/cobUg/ZIRzDxpuoHH3j7iI4X7AUmQXIh6lZxSyIW3lR1IzoAAIu5asgRoPhGiECNhTGAJeMsLmpf2JSooF62WGmMa5NjgRHmkCAxQ1PySPs7NFeBVnukvod1r37YPytXOM1K0QFM7BUSATPnlxPF5ExZt623G53nkBCbbo/2y5xDE2sgHBp8BkPoKO0BxZ0SHYB88GEJPvBXpX4awCTF7TRIDgiO8RAdFV9cY91kIACoalzZUd04omoa9wdZe/OsYZBxSiWMhJNG2JR/qP2M27TCS/frvKzqBRr0wcl9qYieEFP2eSyx5cRPFBGcuz1Pbr7B+m3G++Gknb1yZYbx12LjCc7fyKCZmFFqrj+ZZxRt7z1ViPAOKrgRwZzsyaVvan3ZD0lvXwzwwdt+ZrgjpyNK4GEA5Amu8X9ZYc7KWMmh8eVEfnFF5Y0zwzzBsaiWSv9BI37JcsRFI/fLnU/dR3YdCC9K5HOTdc4fMgSt/Fwaqz9e+iWnD6TmHmm3xIZg4c1SH2QI62aZouxRadln03N3HKf9jIo4R33o2aBVt+7rNvE7H/2pHR1/2qCQXD/YIA4UWISBBzmQXAUiCyYLogaoDKghsLF+LH0pS6LWTyXn9youR1ycUW2jiRDt9yG1Av3yKqLTVmwa3BsLBgWiEwUn5mpkjYnpyhwTHVrCDNt8hwYVZZ+KTX11AwjJGEM+OB4m274+tX3LdySouZ7rfvh3XjXumgUtJkhySwv0iehtRHQfEX1atb2WiP6WiD7pl+9W215NRHcS0WeJ6PKxHRGSk19ysm58fJm/EpZbVWf9hMd+El+ZRo8aN+FvJxUgQ0q5xN6O3w2ZFyQhRb30tYftJr/PUFvfEt1bQthyLRj1XKOFQuAlt81WknfYt098nrB/ct7ovpO27r/JCEbk5DNFRklvV8H93i2F4FYFwf24jln2EmOU3NsB/AKA/5y038DM0VAXRPRkAFcDeAqArwbwO0T0JGZuMIAOuYjaaJow2xFbBlUG7Ieb5qZVc2HqPDkuOXdHxQHRC6GJrkM6sq8+Ft39dduQggm36NWD3jUIvSRIkUNuX63otIIL+1Oi7KQB3CWQ8BipSypJ34b8YEv7vaqJul2Cyz3oHvzezT85sjMFEU5Vnxwzf4SILh55vucCeBczbwH4SyK6E244lD9aeh1lrpKoNHZkJuYqNxZkfNDBkAtCNBZkLLjRRJcqI4pNUp+8OhRk6DNTO+Ys0FWCqq37PGUFLeNw96UP+/QgS5YE90w00cl24TExk+WaPlqht4Xo6pI+9N5jn1Pfd0q7EoRj0/3Im81k2/XcS5Qj2KWpKvLcM/j9DxSC2xEmSHI7GU/uR4joT705e45vuwDA36h97vJtSxERhECXxYThpm08ukMySkb2PEBMRElb7/cc0vOnBJeYYh0zM7Nfn8IcVETJMfq+OFWtqn85/2PaNynLYkrN0GShzBLyDeP2UfeWkFiI7vYQXAd9pKvN5WTfnQYkChJkg1uZZQ+xXZL7ZQBfC+CpAO6BG6ET6P3T7YKIriWiW4no1uahB1fvQc+8n9GFtQ8s7V2GALrn62lHnkBScovOkZJs5ntvP3PE2Idke87s1u2pMl3mV1y65PqNEWSSDi8eiI1bgtNDP2kXRc9z6BAbut/TYwp2hgMzufTACJ13AbhQ7fo4AHf3nONGZj7KzEerM87cTjccJKXEtN9zpNVROuH4/Gm7+XV9plHPdfwxn/n3r4j27+zT6dtwv9Jt2WMz+w2t5wMZO1+6960725rL8kdP6pc+OKcVwS19UUa+PNmyskJw68FBUXJEdL76qkfofB+Aq4noEBE9AcATAXx8+72jtmwmFIMbPweBJipqyW3ZHY34Yx5MfehTfPq4AYId3HeNL1qvchogyE4/lynHMcpS796XBqJ8qW6kZ79uRdH57RJl90tUejYUkMioTa3o5Fl95Lf/7fibKehCfKinWnR1lRE6mfl2IroJwGcA1ABetiyy2r2g/6OTBYArtESH2EjIzbiJU7jykxirIvsh9OWzLUPWvM0cJyrujje4zye/5oZ234wDXVdDDEp67bfqXV/yc6n7oNaz1172TGR73yXTPmZUnDZNUzNVCE73iwmdpO/efgFxMKTn3gvWgAk+yzHR1ZVG6GTmNwB4w6odiRie4FJDqsoV5Hv/G8n3jRlQVYCsb8zcaMJV7BRfFYMKrm//DNlFJqpC2v7kV9+wmn+ih9CyZJFpTz/TbWPGthuD7nm5c13yw6uTle8Mqlv1RrZVbyHCGl2EAkFxSJ6RbarLUZ2yyDiAwG7wBt+Xj7yvqLh1YIopJDuJrq4N4ddaqQk2BmSMIzaZFSqoNAIZA8wqp94q9+kmUNnXW+kluOy+P+v3HdPnHRJcZNItI7gV/SgdX5k2PbUi0+ZKQnCB2CzC1JPUqO9qiVSf7eZHRqSa7XBMfv/tt35iYOeCldD3tzPyb2k3MAmSA+AnTvHrhhxhVZVfrzzZVa5NFk90bjw0SVdYosh28QEH0loBt//ciGNyBJcx9zok1rNPluAyf4CjI2XJQqzIyLb7hGNsouAauEUTnF83jatmgSy2PT8pX136Qxk9N+UHjaPZJdqwVowluD0mucnUrmpzNVQoVAaYzVyFA7MnNm+yzirwrHIqbmZaM1XU3Cn0RyxE940/kRmup09pIXmplbrJ+buic2S+bzda2TEjB0g0Vm2erBrfFhEchzYIYYYLetHmR+yV4bbYMkBwZit7cxR6BGTVRfE7EvDR9xQVty4QpmmuToPkOCY5N9EKOaVWVUBl3d/obAaaVYHgMKuAmQ84VDIIZt5klREudgu3v3F1FdeHZaZkzs+WJbicskNmXZ8r3DugggAAG8tJREFU9GHEX2tOYWb6GKWAWEV0jVq3DKo5IjqEe2pPzhKPIjgSM+KRQ/DPyT80WUSDpkZdXxbcKdgWpvhMp0FyiP/4XSqIAc+camPLAFtPcLNgrvLMtIvM1J5JEei95pqIb5TJuSqG1FenrSW4aAy4LNl0z9s3KrFG7x9vjuii63JEvNoURWSyOtMUkZLL+Nv8DyDLb5mFIzzjhttCRW4G+zDiCrI/egDwB//3+qesPO1RSK4H8qvetGZEUHLGOHJjZ7oGM3VjBt7QpipCtHPtZiohH0nFmgluCZEMRS1T31iO3GKS4/y5h/qB4X37iA3oqrcQfNDBhYUNhOcGaJDABSMM31S5qSthCAzjZher4IZ810Qn5GidZBO1F/4dJ/gyHghM8LlOguSIgWrOMAv3Sw/AmaCHZoDdBNXO5uANT3KbjuDshpuI2lai5hCmJtx2nlzUsfw+n75+F5SbP+83XXdDS1TefI8IDAg+qg55iTnYQ2q6iiAcD0RkFDBgrvb6B8O5uO2/+g4GzMKq4AO3yq22nuREzdkwSEN7YXKRduXO4I0qSh2SbWzcHL6QuSv80PrCb3/47qLi1o7c39EEMAmSA8PlSDXtiyEmK2ZV++OwMQvmqRVfXKLidAlRCGBgHKENfvf41Ft2h+ACkj+SNCVDt3UiqB2Si5Vep1QKQLb+c+APNU9wiSpsYmILPrmGYZqE3ES51Ra0aNwE4n6cwNYv55UckRs92hiQ9TOzkVN2pMgN7CfMrgHMZGjkVtH9YTFTdw+F5PIgAKYGTO1fAkb7a11Vjm+IHMFtVLDy6y1/1HrkCz0wY9/FxrTtJzwxdAguZ3ombXpootQHJyaifNfXkvUceodNAjLqkNuBS3WKh0RTxSQVcrMWVHv1Nl+EWdqCipNh8I1TcCRRd++zNVCuDevyJ7kiMBsYw7CwLdE1TtEV7B7KlIR9YE9wtcp6B3x01YKNmKuK4GYmJrpExQHt9xyJ9Sm7Pt/bp27YZQUXOpDzryUEZ5ElvGAGZtoB5RdLSC2XPNv7x6qtx0TB6T5A98PntKFhmNq25CYKzlo3cfh8EebZZWt9jaqqCiT/o+ej7tQ0YGt9/qQJfy+YGTddJbmJxwPRicO3YNdQzNU+BJKDJzoOioON1K2il+A4o+JClDVK/mzXw6X133yPyvvT/7BHBJdiQMGl6k327frkcqSnzFgVjR1SbPm+ybpSh7oMixGCCEG9LZpEwbnJiqhugK25I7e6dvN7yPD2oeZU0opMmzN5aLP9Xjn3BtXOpWEqmbXNEx0MPvL/lry4XYOyHHYKInobgO8BcB8zf6NvezSA/wrgYria+auY+e+XnWsSJEcMVAuGWTDMwrosd+t+icnPxs6EmNwUwSEit1bFdYgOQN8wS33DKO01wcXVARmCk9KolLgk6z81FdN9UlJL/jD7U0W6G3S1gdtHiJZDcm8gOCE2Ibm6aU3TunHqbT4HFjW4rt1o0E3j0odEzREBpEr9qsoR4mwGqkyUJE51BUME3qh83ytPdAW7ivUpubejO+3CqwB8iJnfSESv8t+XDuU8CZIDM8zcmatmbt2Q5rWNfHNcGfBMT7Ci50qV4ANF5AZAmbCISC0yZZN2JsKf/Px+qjdkAwraV5dGLqMyqR7llg0IAN3gQ25AUo9u0i/H6yEHzpddadPUWtBW3ZKbN1O5rh3JPXwSvHAEd6J+16jHddmh57m5PzY2PMnN2nxKQ4C1MOz+zH/3d8u0gruJdVY89Ey78Fy4EZEA4NcA/B5OFZKjBth4oAY1Fubh2vlofESNN2fgysBuGjSbJvK/BVLzqSNCcFL5wNH32KSN29Gehwh/8gv7RHBQak2bnjYmNtO0Ck+PvUaKyKhBbJpqAoyIUJmuqVobEj5RpJcVyblZ01JyQ9242dS25k6hCbGJclvUo4lN45atd2TbL9u8BuaRj4A5fBg4tIlb7nzzyucuWB008AO5BjyWme8BAGa+h4geM+agaZAcM8yWJ7kt73wGvOlhAVTQ8wekBBf54VLlluwfbVPrTMAnf/GVe3XLvfiTn38Fnvqy41n/Wy5VJDtkUarYAgFyRGyB1HJR1xRDEVYhN0YcLfXkps1S3tpaC7ktwy3zd+KyzWtA8wU++I9vW/v5CzJYzSd3LhHdqr7fyMw3rr1PmAjJwXqSqy0wX7gXgshNP7jhiw8zhBUptIjolKrT+0fnAD5x4/6TWhY5M1WTmI2/t9tVwTur7cFfpxUXlN+Ooz/O0RM5R0nJ3v/WuMTeSL01nuDmC/DWHFgsYOeLXSE3jVvm79zV8xd0sYK5ej8zH13x9PcS0flexZ0P4L4xB02C5IgBOrkALWpnzlin5AgA2w0VXcsTHBKCy03GIt9vfdtEiS1Fot66hJZ8F/JrWuJq10XhcZbYyCvnwaqHXKlc8MdxMLFlsm80tiW4+QJcN8BiDj65BV7MccK+e62Pq2Ai2FVrFe8D8EIAb/SfvzXmoEmQHKwFndxyv/YnT7oXQ8aQa2KVkSM4TglO+do+/p+v25972gGI4yUivD6CE19cappKjpqKdkb+M6XCArK+OL89HfgglJ5xIDuX32ZVYMET3HzR60MrOBhYV+ChZ9qFNwK4iYheDOCvAVw55lwTITkGP3TSR9gedmaqj5TB2m6yatYEbcu7PvaOU4/YIqTE5kmpJb08wUEpubYgvlVv1NgssUm98FC9atu3YUUnZBcI1VrAupw3XtQ7fzYF08aaSK5n2gUA+M5VzzUJkuO6RnP//TjR3BS1X/HoH/Z+nXakWcAR3B+85+DWHw764SxgGkVuyYCTaW0oPLmFoICM8KGCDpQquT4zNZqmsUt2nfNoNA24Xmz/oRRMH+JDnhgmQXJPuuQJOHHrTd0NNlEIHgeZ4AC0gQe/Hi3otkfmbY7gJCm3trH/TJ5po/4yU5LSZKamhyQ1/FHvsfocPpG34OBinXly68QkSG4U9GCQpwO45w8mNd3F9BSlq6shNMH5ORICuckxQAj0hPNrRBNCtyQHQ0Af0aVlWKFKoWeY3oKDg7GR+T3E9EkueYlOh5mVNLmFyVqQtIkaU8doAgvHsY6qio8sNllphJJjTWrMaMch53g+3DQZ1I8Bx1KKVXCgUZTcqpBpCGWG8//v4M+NeclLjufHCVAmaa4dQKv+GC2RSaJvkyE4HSDo+wX2xEYNB1ILhCeDJzSK/Dod9/9+lQE2NrbzSApOFaQulYlg2iRHJlYJBxyXXHs8rI/6RUzNVqCr4pSyg4p4UkJ2API1q6mCU4THsK3S7gtchBsyRcmdBphi4GHaf3Uy1DVh6j3dMS659vj6fwVFxclnTr0J+elBKvXi007Cdmuj81DmuHB+gZHJwQ2Oza5e800WTAk6C2Bo2UtMW8lVVXByf+jDr9nv3uwKtHoDWgUns1G5fEBGp6zNOL9cmLVK5QsGtjSIE3tzJCaJuwA+8Df/Z28/n3Phj/mZsVr3gbgTSNr1dfT1BJWB2dzAseqqTrpQwQGABMEmhkmTHBFN0cTfVej5QPUkyBGJ+R2YCEQcEqIpnMNNvAwgJr0ev9kQufXt85yLXg4/DmXmJhLzV6K33o9XoqwHFyXwsCLYvywn/vCn97knu4NUxWXhOSodXIBS4jMA23ab+5ToJxBVKhDhA39x/Y76/oG//g8AgOc8/hWt304jZxIXHHxM8J950iQHy7j5U6/f717sCrIEpwhNTFAGVPItOiZrkHuB+PzBfnJl9pO/COnd/OdvWut9fOCvbgjrz3m8GocvkFveAXPMXFmK9A8YSjLwdqAnMTlAuOQlx5HNEwE6yk3MVTbk1JDU6+p5WbxS4wpu5ir2uWuW3fDxbkYX3HzHz+7qfUWEd9HLwzozO5lp+0mv4ABAgk4Tw6RJ7uYv/sp+d2HtuOQliYJLyU6IK5AbXBWD7Kv8cmwUIxpuj/UzyHNlQGzdrPN7HJ4Wc7bgNMP0OG7aJHeQ8LQf9uSWG8RDt0nllAQmhezYmZuu/JNBltqQPLlPJ97cdHxErtLBEkCWcMsfH0y/ZsG0UMzV0xRPe3GPeUp5ggNaUabXRdVxJMrE6U+wM7c91mxpHklBwS5BKm0mhqU2DBFdSEQfJqI7iOh2Ivox3/5oIjpBRJ/zn+eoY15NRHcS0WeJ6PLdvIFTBvrfPiU8atv0qMZhAFADP1kPwqQ90ffK++OMi7JaPW1jRfidj5ZZqgr2CDxy2UOMcdTUAK5j5m8A8G0AXkZET0Y7B+ITAXzIf4ffdjWApwC4AsAvEVFJjEpB8SdrVacJjloy0+tsCLYi2EoRXwXYmSM2WxHsjPC7v/vqvb+3gtMW6ajWfcteYinJMfM9zPwJv/4AgDsAXAA3B+Kv+d1+DcD3+vXnAngXM28x818CuBPAM9bd8YOAaPLrJD0kTLFolLJTqk6Tm63gF6/eZggqrqBgL+EGeF2+7CVW8sn5yV6/BcDH0D8H4gUA/rs67C7fdlri6T90HGQQ+94U4ukRJeqgvgeHHMdSX9xsSRWErVyAwoDx39578EdtKZgQ9sEUHYPRJEdEZwF4D4CXM/OXKTd7k98109a5dSK6FsC1AHDRRReN7capC0KnagEpwSUmbFhnRXhhQEq/2HhH9pnEH333AR89uWBycMnA02O5UclTRLQBR3DvYOb3+uZ7/dyHSOZAvAvAherwxwG4Oz0nM9/IzEeZ+eh555233f5PGs/4wbdk26PAgpimneACWjN05haeOZ+b3SD/CTSbeiG/7O19FhQE2JHLHmJMdJUAvBXAHcysM1llDkQgngPxfQCuJqJDRPQEAE8E8PH1dfnUwLc+3xNc+sPWE1lNVZ4EGiJfXRRR9X43I5UObhHfXEHBfoDCeIbDy15ijLl6KYAXAPgUEX3St70GPXMgMvPtRHQTgM/ARWZfxswHsz5rGXSyG+L12BeXWU9z6DRkvovgr5PaVHYVXiUtrmA/cKr65Jj5o+ivtMzOgcjMbwDwhh3065TGt13zFqeRh+pTkfjhkPfZZYlOyr5ksXIcwYLxyV955XpupKBgJZTa1dMG8az3DFZMpnPhwveeoERU2aDJME2sJPXZy6wFBXuACQYeCsmtGc+86vql81JERNdHcEuUHPmMEvgaV0pN44KCvQZjknM8FJJbNyTNg2nQPxHNQD9EcD1Ex2gVoxTzMx34qTAKpo6i5A42nvV9bwb5SgUopRVA3SUq10rbUpJLyE7Or6ckvP1Nr0BBwb5hehxXSG6dcBO6+IRcKJOyB1FZF5CUdsWfSPeVQzSZTvAPrOD0Ak1wUNRCcuuEFCMw/PwUI5xkGVLL1bTqNJNOEEJdu6Bg3yCR/jWBiL4A4AEADYCamY9u5zyF5NaEf/acnwM2TEw2q4y4kBJYatL6Nkd26qRK2k1xwMKC0weEXUn0/Q5mvn8nJygktya4IWR8uggnpmpObfUot1TFRf45f6F4oE1NeGu9pYKC1THBwEMJxq0B3/Fdbwwjom47hJ6aox2/HIMNh+GWcguXf82C/UZuAvPcMvJsAG4hotv8gB7bQlFyO8R3/rN/72bECqOEeEXH1M7TAHTV3GB6CRzR+YgrG87m1IX9468FBfuD1Xxy5xLRrer7jcx8Y7LPpcx8tx/G7QQR/Rkzf2TVbhWS2wG+61mvB4yQWT4vjrTpugLRAYn/bUlKyfSMhILTEStEV+9fFkhg5rv9531E9Jtwg++uTHLFwNkBwiinVkZfcJ8h4CAyTgUgaBnBJUSWTSXpybcr5mrB/mKkqTrCXCWiM4nobFkHcBmAT2+nV0XJbROXPf21IGP8YJZS4eA+yc+olZqr8l2UXaTyclAqTsDquw5GFBTsO/QP+87xWAC/6QfnnQH4DWa+eTsnKiS3TVDDYFhQA4AM0HA0V6r2y5FlgCgqrJftwDKiU5+a7MSULc64gilhTXlyzPx5AN+8jnMVktsGrvjGfwcYA4KXa9ZFPt2kznDE05BfZxBTW0Tvl1WSeIPpCsS5cmldayG8gn3GFIc/LyS3HdQNaAY3OGVDgFEzETF11By0miNEpmvAGF9dWE8ILlF5BQX7hkJypz6e8/hXgCoDthbEsxB0YAbYEAwDmBGYTSA7mQDagkEyUCYBZMSX50++LBckV+aVJAoXFOwbmIGm1K6e+qhrADMQEZgaN+CIAVATaMP55gwBFhaAcbMFqkCE9tu1UdgeVZdBdiw6oPjnCqaBouROfXDdRLEAJoJMz0i15HAY12YZ1ADEBLIMNu6TDLXjwW2nDymZFYIrmAomSHIls2oFXHbk+cBiDq4bp+jqBrTwn7UFNdZ9WoZprAtEBF9dHHhIl7DdY1sESMDFv3z9Gu60oGAbYLjyxjHLHqIouRXA8zm4qpyC41mIJLnvDJr53wzjarKIXLSVGpdlQtZ/KsLrNV/34wYLCnYExhSniitKbgWcaG5yRHdyC5gvwIsFeD4HFjVoUTsV5xdnqrIr2A8kxr0qrgNGxi7N7JsxVR//H9+Mx//HN+/4fgsKVgLDBR7GLHuIouRWxC3zd4b1yzavAW1uApsbwMam88M1GwC7+U8tAcYA1Bgn7kJk1aWTWILy1SEOJgj5+fbg6iCv/oTYws7ytWXBx/9KS3R/9dKf2J0HUlCgUXxyBwu3zN8Je3ILPF84H521gHW+ODS6rhWDai5We7E5C6h2JG0rQBNeQcGuYb1DLa0FheR2iBP1u5wJW9dA44IPsBZkk8CDZVXE7wv6beybywcm2kE4AwaIjrJzGDoUM7ZgdzGS4PaY5Iq5ugZwXYMWC6BpAFs5n4P8Y0o0iSkOQFi1bhiwFHLpJPAQlYEBEamJyerqXqndOCKVRCKwX/jff3ydj6HgdAfDWTMTQ1Fya8AJ+25wY8HWBoILJqsoNqvN1J4AhDJbY1WnBuRE/H1QzfX8YIrf7uJfuh4X/2JJOSlYI4qSO7jgpnFKjhnkCY+sBRrjzVSfCGwBsqSIrx2aKRCdlGwpwuP0k3rUHMMHJ8iRmQpYZFGSiAvWBi5lXQcatvHhcW+y2vZXi3w6CSoKfriwuOovR3gmHY4JgPVEZhVpibpDD4kxuuQ18ON58S9cDzDhCz963ZoeRsFpCQZ4gnlyheTWhBP23bj8kS8CFgZUVSBTA5UBjAXVDCKGadiljRifTUIE0/hUEwCWGAYUhuSiKr2KJzYDn1oSp4+ACKy+E9NKSu0JP/+WyGz+/CsL6RWsiD2uZhiD4pNbJ7y52n5an1bifXM+OZgaRH46IZXw2RNpDcovSktJ/HP6O9Cv4JZtB/C11x/f2fMoOP0wQZ9cIbk1QvvldM4csQtCSNCBWKWTdIiOWyLTpm0ISlBEgJS0yfe+GtksAeq0k+Tv7+veVIiuYCTU3/3SZQ9RzNV1omnAvmyFVBoJNQyqPIE1AFcUiE0IjA3iAEROyenqB3bWKvtUFGapoRWjVpmyfX479RklHyf4up87DrLA5179yjU+rIIDiVLxcLARlJxto6whZ04nA3uTtSU6jkzWNgobEyF8IELaYdsqCLKxeluq6IBeUnPbEnXHwJNeX1RdwRAY3DSjlr1EIbk14kRzU0t0TeNkuTdVxWSNKiCUqUrqM82VI6/u2+0Ubw/7dU1ZR34qb65DdDnzt/8ev/5nbsDXv/aGNT+5ggMBxiSHWlpKckR0IRF9mIjuIKLbiejHfPtriehvieiTfvludcyriehOIvosEV2+mzcwOXiTlRvr/0FtG3QQAtHkFtY5pIyEdgukQQkJOkQEmJq1me+B7JIl2lfQR3yKXP+Xny5EV5AB23HLHmKMT64GcB0zf8JP9nobEZ3w225g5ihlnoieDOBqAE8B8NUAfoeInsTMe6tR9wmhxGs28+kkBlg0fgRh50czZADj0klM49I8JCNEskMAuPSScGJ1kRoh546NI0kYeNu1nUNCD8PE/tpxZ91HPBgAtWQLZRaHBGa05FtQoMAA+FRMIWHme5j5E379AQB3ALhg4JDnAngXM28x818CuBPAM9bR2VMBUuKFpgFLGgmzK9gXxZaMNRdGDhZfXKLUInM2VXZMcXukBskl4enzpL4+1uehjq8vl9JS5sspyIJ5kkpuJZ8cEV0M4FsAfMw3/QgR/SkRvY2IzvFtFwD4G3XYXciQIhFdS0S3EtGtX/ziF1fu+KTB1hOcypfz4XUxQ6EIDzkCs5r8MkEIHYBIiQ7oMVXRXSx1fXkAYmIkdS7Vx4KCBFMMPBCPDPkS0VkAfh/AG5j5vUT0WAD3w70arwNwPjO/iIh+EcAfMfOv++PeCuD9zPyegXN/EcCD/nynAs5F6etu4VTq70Hr6+OZ+bztXoCIbvbXGYP7mfmK7V5rFYzKkyOiDQDvAfAOZn4vADDzvWr7rwL4bf/1LgAXqsMfB+DuofMz83lEdCszH12h7/uG0tfdw6nU39LXGHtFWqtiTHSVALwVwB3MfFy1n692+z4An/br7wNwNREdIqInAHgigI+vr8sFBQUF4zFGyV0K4AUAPkVEn/RtrwFwDRE9Fc5c/QKAlwIAM99ORDcB+AxcHPBlp0tktaCgYHpYSnLM/FEgO5bF+weOeQOAN6zYlxtX3H8/Ufq6eziV+lv6egpgdOChoKCg4FREKesqKCg40Nh3kiOiK3z5151E9Kr97k8ORPQFIvqUL1+71bc9mohOENHn/Oc5y86zS317GxHdR0SfVm29fdvPkruevk6yPHCgnHGqz7aUX/aBmfdtAVAB+AsAX/P/t3f+rlFEURg9FyEWwUZFSRnFJpVsYaOkFJJG/4MUgo2NhUUgja2CtYVVsNBGxJSCja0IJmrhT6wMSWFhqci1eC+Im3mrBp17Hb4DwwyPKQ4fu3d3mP1mgSlgHZiLdGp4fgAOj61dB5br8TJwLchtHhgBL37lBszVjPcDszX7fcGuV4ErHedGu84Ao3p8AHhdnbJm2/JNmW+fW/Q3uVPAW3d/7+5fgLuUWtj/wDlgtR6vAucjJNz9MfBpbLnlFlq5a7i2iHZt1RmzZqv6ZYPoIfdbFbAEOPDQzJ6a2cW6dtTdN6G8wIAjYXa7abllzXvP9cA+GKszps/2b9Yvh0D0kOv6aUrG272n3X0ELACXzGw+WmiPZMz7JnAcOAlsAjfqegrXWme8B1x298+TTu1Yy+CbOt8+iB5yf1wBi8DdP9b9NnCf8rV+a6f1UffbcYa7aLmly9vdt9z9m5f/srvFj0umcNeuOiOJs23VL7Pm2xfRQ+4JcMLMZs1sivIcurVgp58ws+n6HD3MbBo4S6mwrQFL9bQl4EGMYSctt3SVu6z1wFadkaTZqn45geg7H8Ai5U7QO2Al2qfD7xjlLtQ68HLHETgEPALe1P3BIL87lMuQr5RP5wuT3ICVmvUrYCGB623gObBBeePNJHE9Q7l82wCe1W0xcbYt35T59rmp8SCEGDTRl6tCCPFP0ZATQgwaDTkhxKDRkBNCDBoNOSHEoNGQE0IMGg05IcSg0ZATQgya76CedvuWGuLMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now lets subset the SNMR data into the Keep River area where we have a water\n",
    "# table map to differentiate saturated from unsaturated\n",
    "\n",
    "wt_infile = r\"C:\\Users\\PCUser\\Desktop\\EK_data\\Interp\\WaterTable\\KR_wt_depth_clipped.tif\"\n",
    "\n",
    "wt_src = rasterio.open(wt_infile)\n",
    "wt_arr = wt_src.read()[0]\n",
    "\n",
    "# Romve null values\n",
    "\n",
    "null = -3.4028235e+38\n",
    "\n",
    "wt_arr[np.isclose(wt_arr, null)] = np.nan\n",
    "\n",
    "plt.imshow(wt_arr)\n",
    "plt.colorbar()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we keep only points that fall within this water table map\n",
    "# TODO get a better solutions - get polygon of water table map\n",
    "\n",
    "wt_bounds = wt_src.bounds\n",
    "\n",
    "df_acquisitions['wt_depth'] = np.nan\n",
    "\n",
    "# Iterate through the sites and add the raster value to the dataframe\n",
    "# if it is within the extent of the raster\n",
    "for index, row in df_acquisitions.iterrows():\n",
    "    # get the coords\n",
    "    x, y = row.mid_X, row.mid_Y\n",
    "    # Check if point is within bounds - HACK!\n",
    "    if point_within_bounds(x,y, wt_bounds):\n",
    "        \n",
    "        wt = next(wt_src.sample(np.array([[x,y]])))[0]\n",
    "        \n",
    "\n",
    "        if not np.isclose(wt,-3.4028235e+38):\n",
    "            df_acquisitions.at[index, 'wt_depth'] = wt\n",
    "\n",
    "# Drop null columns\n",
    "\n",
    "df_acquisitions_Keep = df_acquisitions.dropna(how = 'any', subset = ['wt_depth'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now get the SNMR inversions from this area, clip it to the water table\n",
    "\n",
    "df_Keep = df_inversions[df_inversions['acquisition_id'].isin(df_acquisitions_Keep.index)]\n",
    "\n",
    "df_Keep = df_Keep.join(df_acquisitions, on = 'acquisition_id')\n",
    "\n",
    "# Create mask for if water table depth is more than depth to\n",
    "\n",
    "mask = df_Keep['Depth_from'] < df_Keep['wt_depth']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a saturated or unsaturated label\n",
    "\n",
    "df_Keep['label'] = 'unsaturated'\n",
    "\n",
    "df_Keep.at[df_Keep[mask].index, 'label'] = 'saturated' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate K\n",
    "\n",
    "df_Keep['K'] = SDR_K(df_Keep, C = best_C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find log K\n",
    "\n",
    "logK = np.log10(df_Keep[df_Keep['label'] == 'saturated']['K'].values)\n",
    "\n",
    "# Romve infinities which are divide by zero errors\n",
    "logK[logK==-np.inf]=np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAGDCAYAAADEegxVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deZhkZX328e8tKKigQGgQAdMYcQHjOiJKNCgSMRhAExWDOhgNWVATNeoQcon6ypsxmrjELRNEUCcgGgkgKpsiMa+CA8omKLyCMIJM476CA7/8cc4ca5rqZZqpqu7p7+e65uo6z9l+VdNddz3nnHpOqgpJkgDuMeoCJEnzh6EgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCloQklyf5Bkj2G8lecgU8w5Lcvawa5IGyVDQopLkhCRv7ZneM8nNSV67oduqqpVV9Qcbuk9pPjMUtGgleQzwBeDYqvrnUdczF0k2H3UN2rQYClpwkmyR5F1Jbmr/vSvJFj3zX99++r8pycv7HQJKshdwLvD3VfXeGXb5jCTXJPlhkvclSbuNw5N8qX2cJO9MsibJj5NcluSRSY4ADgNen+RnSc5ol39EkvOT/CjJlUkO6qntt5KckeQnSb6a5K3r9tPOryRHJrkGuKZte3eSG9t1Lk7ylJ7l35TkE0k+luSnSS5P8tAkR7X13phkxh6PFgdDQQvR0cDewGOARwN7Af8AkOQA4DXAM4CHAL/fZ/29gM8Br66q42axv2cDT2j39XzgmX2W+QPgqcBDgW2AFwDfr6oVwErgn6pqq6r6oyT3BM4AzgZ2AF4JrEzysHZb7wN+DjwAWNr+m+wQ4InAHu30V9vXYzvgP4BPJNmyZ/k/Aj4KbAt8DTiL5u9/Z+AtwL/N4nXQImAoaCE6DHhLVa2pqgngzcCL23nPBz5cVVdW1S/aeZPtDfwY+Ows97e8qn5UVTfQHG56TJ9lfg1sDTwcSFVdVVU3T7G9vYGt2u3eXlWfBz4NvDDJZsAfA8dU1S+q6hvAiX228Y9V9YOq+iVAVX2sqr5fVWvbQ2FbAA/rWf6/q+qsqloLfAIYa/f/a+BkYDzJNrN8PbQJMxS0ED0Q+E7P9HfatnXzbuyZ1/t4nffRfLI+J8m2s9jf93oe/4LmDX097Rv7e9tt35JkRZL7TVP/jVV156TnsDPNm/Xms3gO67UleW2Sq9pDVz8C7g9s37PILT2PfwncWlV39EzT73lp8TEUtBDdBPx2z/SD2jaAm4Fdeubt2mf9O2h6GzcAZ03z5r1Bquo9VfV4YE+aw0ivWzdr0qI3Absm6f37exDwXWACWMvMz6HbZnv+4A00vaRtq2obmp5Q5v5stFgZClqITgL+IclYku2BNwIfa+edAry0PZF7n3beXbSHTZ4H3Ap8Jsl9705BSZ6Q5Int+YKfA7+iCR9oPqU/uGfxC9tlXp/knkn2pTnmf3L76f1TwJuS3CfJw4GXzLD7rWmCZALYPMkbgY0SdFp8DAUtRG8FVgGXAZcDl7RtVNVngffQHPu/Fvhyu85tkzdSVbcDz6V5Az8jyb3vRk33A/4d+CHNoaDvA+9o530I2KO90ui/2v0eBDyLJpTeD7ykqq5ul38FzeGf79GcHD6pX/09zqI5P/Ktdt+/ov8hJ2lG8SY72pQleQRwBbBFe5J1wUnyNuABVdXvKiRpo7KnoE1OkuckuVd7EvltwBkLKRCSPDzJo9rvPuwFvAw4ddR1aXEwFLQp+gua4+v/n+a4/l+NtpwNtjXNeYWf05wj+WfgtJFWpEXDw0eSpI49BUlSx1CQJHUW9AiL22+/fY2Pj4+6DElaUC6++OJbq2qs37wFHQrj4+OsWrVq1GVI0oKS5DtTzfPwkSSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSpYyhIkjqGgiSps6BHSZU0XOPLzpx2/vXLDxxSJRqUgfUUkhyfZE2SKya1vzLJN5NcmeSfetqPSnJtO++Zg6pLkjS1QfYUTgDeC3xkXUOSpwEHA4+qqtuS7NC27wEcCuwJPBA4N8lDq+qOAdYnSZpkYD2FqroA+MGk5r8CllfVbe0ya9r2g4GTq+q2qroOuBbYa1C1SZL6G/aJ5ocCT0lyYZIvJnlC274zcGPPcqvbtrtIckSSVUlWTUxMDLhcSVpchh0KmwPbAnsDrwNOSRIgfZatfhuoqhVVtaSqloyN9b3FqCRpjoYdCquBT1XjIuBOYPu2fdee5XYBbhpybZK06A07FP4LeDpAkocC9wJuBU4HDk2yRZLdgN2Bi4ZcmyQtegO7+ijJScC+wPZJVgPHAMcDx7eXqd4OLK2qAq5McgrwDWAtcKRXHknS8A0sFKrqhVPMetEUyx8LHDuoeiRJM3OYC0lSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHUMBUlSx9txSurMdLtNbfrsKUiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOt55TVpkvLuapmNPQZLUGVgoJDk+yZokV/SZ93dJKsn2PW1HJbk2yTeTPHNQdUmSpjbInsIJwAGTG5PsCuwP3NDTtgdwKLBnu877k2w2wNokSX0MLBSq6gLgB31mvRN4PVA9bQcDJ1fVbVV1HXAtsNegapMk9TfUcwpJDgK+W1WXTpq1M3Bjz/Tqtq3fNo5IsirJqomJiQFVKkmL09BCIcl9gKOBN/ab3aet+rRRVSuqaklVLRkbG9uYJUrSojfMS1J/B9gNuDQJwC7AJUn2oukZ7Nqz7C7ATUOsTZLEEHsKVXV5Ve1QVeNVNU4TBI+rqu8BpwOHJtkiyW7A7sBFw6pNktQY5CWpJwFfBh6WZHWSl021bFVdCZwCfAP4HHBkVd0xqNokSf0N7PBRVb1whvnjk6aPBY4dVD2SpJn5jWZJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUmfzURcgadMxvuzMKeddv/zAIVaiubKnIEnqGAqSpI6hIEnqGAqSpI6hIEnqGAqSpI6hIEnqDCwUkhyfZE2SK3ra3p7k6iSXJTk1yTY9845Kcm2SbyZ55qDqkiRNbZA9hROAAya1nQM8sqoeBXwLOAogyR7AocCe7TrvT7LZAGuTJPUxsFCoqguAH0xqO7uq1raTXwF2aR8fDJxcVbdV1XXAtcBeg6pNktTfKM8p/Bnw2fbxzsCNPfNWt22SpCEaSSgkORpYC6xc19RnsZpi3SOSrEqyamJiYlAlStKiNPRQSLIUeDZwWFWte+NfDezas9guwE391q+qFVW1pKqWjI2NDbZYSVpkhhoKSQ4A3gAcVFW/6Jl1OnBoki2S7AbsDlw0zNokSQMcOjvJScC+wPZJVgPH0FxttAVwThKAr1TVX1bVlUlOAb5Bc1jpyKq6Y1C1SZL6G1goVNUL+zR/aJrljwWOHVQ9kqSZ+Y1mSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQwFSVLHUJAkdQZ2kx1JozG+7MxRl6AFzJ6CJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKkzsFBIcnySNUmu6GnbLsk5Sa5pf27bM++oJNcm+WaSZw6qLknS1AbZUzgBOGBS2zLgvKraHTivnSbJHsChwJ7tOu9PstkAa5Mk9TGwUKiqC4AfTGo+GDixfXwicEhP+8lVdVtVXQdcC+w1qNokSf0N+5zCjlV1M0D7c4e2fWfgxp7lVrdtd5HkiCSrkqyamJgYaLGStNjMlxPN6dNW/RasqhVVtaSqloyNjQ24LElaXIYdCrck2Qmg/bmmbV8N7Nqz3C7ATUOuTZIWvWGHwunA0vbxUuC0nvZDk2yRZDdgd+CiIdcmSYve5oPacJKTgH2B7ZOsBo4BlgOnJHkZcAPwPICqujLJKcA3gLXAkVV1x6BqkyT1N7BQqKoXTjFrvymWPxY4dlD1SJJmNl9ONEuS5gFDQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUMRQkSR1DQZLUmVUoJNlnNm2SpIVttmMf/SvwuFm0SRqC8WVnjroEbaKmDYUkTwKeDIwleU3PrPsB3kNZkjYxM/UU7gVs1S63dU/7T4A/GVRRkqTRmDYUquqLwBeTnFBV3xlSTZKkEZntOYUtkqwAxnvXqaqnD6IoSdJozDYUPgF8EDgO8I5okrSJmm0orK2qDwy0EknSyM32y2tnJPnrJDsl2W7dv4FWJkkautn2FJa2P1/X01bAgzduOZKkUZpVKFTVboMuRJI0erMKhSQv6ddeVR/ZuOVIkkZptoePntDzeEtgP+ASwFCQpE3IbA8fvbJ3Osn9gY8OpCJJ0sjMdejsXwC7b8xCJEmjN9tzCmfQXG0EzUB4jwBOGVRRkqTRmO05hXf0PF4LfKeqVg+gHknSCM3q8FE7MN7VNCOlbgvcPsiiJEmjMds7rz0fuAh4HvB84MIkDp0tSZuY2R4+Ohp4QlWtAUgyBpwLfHIuO03yauDlNOcpLgdeCtwH+DjNSKzXA8+vqh/OZfuSpLmZ7dVH91gXCK3vb8C660myM/AqYElVPZLmxPWhwDLgvKraHTivnZYkDdFs39g/l+SsJIcnORw4E/jM3djv5sC9k2xO00O4CTgYOLGdfyJwyN3YviRpDma6R/NDgB2r6nVJngv8HhDgy8DKueywqr6b5B3ADcAvgbOr6uwkO1bVze0yNyfZYYqajgCOAHjQgx40lxIkSVOYqafwLuCnAFX1qap6TVW9mqaX8K657DDJtjS9gt2ABwL3TfKi2a5fVSuqaklVLRkbG5tLCZKkKcwUCuNVddnkxqpaRXNCeC6eAVxXVRNV9WvgU8CTgVuS7ATQ/lwzzTYkSQMwUyhsOc28e89xnzcAeye5T5LQDK53FXA6v7lvw1LgtDluX5I0RzOFwleT/PnkxiQvAy6eyw6r6kKaS1kvobkc9R7ACmA5sH+Sa4D922lJ0hDN9D2FvwVOTXIYvwmBJcC9gOfMdadVdQxwzKTm22h6DZKkEZk2FKrqFuDJSZ4GPLJtPrOqPj/wyiRJQzfb+yl8AfjCgGuRJI3YXO+nIEnaBBkKkqSOoSBJ6hgKkqTObIfOlqS7ZXzZmdPOv375gUOqRNOxpyBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqSOoSBJ6hgKkqTOSEIhyTZJPpnk6iRXJXlSku2SnJPkmvbntqOoTZIWs1H1FN4NfK6qHg48GrgKWAacV1W7A+e105KkIdp82DtMcj/gqcDhAFV1O3B7koOBfdvFTgTOB94w7Pqk+WB82ZmjLkGL1Ch6Cg8GJoAPJ/lakuOS3BfYsapuBmh/7jCC2iRpURtFKGwOPA74QFU9Fvg5G3CoKMkRSVYlWTUxMTGoGiVpURpFKKwGVlfVhe30J2lC4pYkOwG0P9f0W7mqVlTVkqpaMjY2NpSCJWmxGHooVNX3gBuTPKxt2g/4BnA6sLRtWwqcNuzaJGmxG/qJ5tYrgZVJ7gV8G3gpTUCdkuRlwA3A80ZUmyQtWiMJhar6OrCkz6z9hl2LJOk3/EazJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOqO6R7MkrWd82ZlTzrt++YFDrGRxs6cgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSeqMLBSSbJbka0k+3U5vl+ScJNe0P7cdVW2StFiNsqfwN8BVPdPLgPOqanfgvHZakjREIwmFJLsABwLH9TQfDJzYPj4ROGTYdUnSYjeqnsK7gNcDd/a07VhVNwO0P3fot2KSI5KsSrJqYmJi8JVK0iIy9FBI8mxgTVVdPJf1q2pFVS2pqiVjY2MbuTpJWtxGcZOdfYCDkvwhsCVwvyQfA25JslNV3ZxkJ2DNCGqTpEVt6D2FqjqqqnapqnHgUODzVfUi4HRgabvYUuC0YdcmSYvdfPqewnJg/yTXAPu305KkIRrpPZqr6nzg/Pbx94H9RlmPJC1286mnIEkaMUNBktQZ6eEjSZqN8WVnTjv/+uUHDqmSTZ89BUlSx1CQJHUMBUlSx1CQJHUMBUlSx1CQJHW8JFUagZkusZRGxZ6CJKljKEiSOoaCJKljKEiSOoaCJKljKEiSOoaCJKnj9xSkAfG7CFqI7ClIkjqGgiSpYyhIkjqeU5CmMd15AW8BqU2RPQVJUsdQkCR1DAVJUsdQkCR1DAVJUsdQkCR1DAVJUmfo31NIsivwEeABwJ3Aiqp6d5LtgI8D48D1wPOr6ofDrk+aLcc20qZoFD2FtcBrq+oRwN7AkUn2AJYB51XV7sB57bQkaYiGHgpVdXNVXdI+/ilwFbAzcDBwYrvYicAhw65Nkha7kZ5TSDIOPBa4ENixqm6GJjiAHaZY54gkq5KsmpiYGFapkrQojCwUkmwF/Cfwt1X1k9muV1UrqmpJVS0ZGxsbXIGStAiNJBSS3JMmEFZW1afa5luS7NTO3wlYM4raJGkxG3ooJAnwIeCqqvqXnlmnA0vbx0uB04ZdmyQtdqMYOnsf4MXA5Um+3rb9PbAcOCXJy4AbgOeNoDZJWtSGHgpV9SUgU8zeb5i1SJLW5zeaJUkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1DEUJEkdQ0GS1BnF2EeStFHNdGvU65cfOKRKFj57CpKkjj0FbdL8BCltGHsKkqSOoSBJ6hgKkqSOoSBJ6niiWdImb7oLDrzYYH32FCRJHUNBktQxFCRJHc8paOT8gpk0f9hTkCR17CloUZupl6JNnz3V9dlTkCR1DAVJUsdQkCR1PKegec/j/pqvNsXzEfYUJEmdeddTSHIA8G5gM+C4qlo+4pI0C5viJyZpMZpXPYUkmwHvA54F7AG8MMkeo61KkhaPVNWoa+gkeRLwpqp6Zjt9FEBV/WO/5ZcsWVKrVq2a8/4GNXLiKD81353n5LF7aeG4O+8jSS6uqiX95s2rngKwM3Bjz/Tqtk2SNATz7ZxC+rSt15VJcgRwRDv5syTfHEghb2N74NYBbXtjmnWdG3m/G2pgr+dGtBBqBOvcmBZCjdCnzrv59/zbU82Yb6GwGti1Z3oX4KbeBapqBbBi0IUkWTVV92o+sc6NZyHUCNa5MS2EGmG4dc63w0dfBXZPsluSewGHAqePuCZJWjTmVU+hqtYmeQVwFs0lqcdX1ZUjLkuSFo15FQoAVfUZ4DOjroMhHKLaSKxz41kINYJ1bkwLoUYYYp3z6pJUSdJozbdzCpKkETIUZiHJ3yWpJNuPupZ+kvyfJJcl+XqSs5M8cNQ1TZbk7Umubus8Nck2o66pnyTPS3JlkjuTzKurUpIckOSbSa5NsmzU9fST5Pgka5JcMepappNk1yRfSHJV+//9N6OuabIkWya5KMmlbY1vHsZ+DYUZJNkV2B+4YdS1TOPtVfWoqnoM8GngjaMuqI9zgEdW1aOAbwFHjbieqVwBPBe4YNSF9FpAQ8CcABww6iJmYS3w2qp6BLA3cOQ8fD1vA55eVY8GHgMckGTvQe/UUJjZO4HXM+lLdPNJVf2kZ/K+zMNaq+rsqlrbTn6F5jso805VXVVVA/lC5N20F3BtVX27qm4HTgYOHnFNd1FVFwA/GHUdM6mqm6vqkvbxT4GrmGejJ1TjZ+3kPdt/A//bNhSmkeQg4LtVdemoa5lJkmOT3AgcxvzsKfT6M+Czoy5igXEImAFJMg48FrhwtJXcVZLNknwdWAOcU1UDr3HeXZI6bEnOBR7QZ9bRwN8DfzDcivqbrs6qOq2qjgaObgcRfAVwzFALZOYa22WOpum6rxxmbb1mU+c8NOMQMNpwSbYC/hP420k97nmhqu4AHtOegzs1ySOraqDnaxZ9KFTVM/q1J/ldYDfg0iTQHO64JMleVfW9IZYITF1nH/8BnMkIQmGmGpMsBZ4N7FcjvBZ6A17L+WTGIWC0YZLckyYQVlbVp0Zdz3Sq6kdJzqc5XzPQUPDw0RSq6vKq2qGqxqtqnOaP8nGjCISZJNm9Z/Ig4OpR1TKV9uZJbwAOqqpfjLqeBcghYDaiNJ/0PgRcVVX/Mup6+kkytu4qvST3Bp7BEP62DYVNw/IkVyS5jOZw17y7vA54L7A1cE576ewHR11QP0mek2Q18CTgzCRnjbomaIaAoTkseBbNSdFT5uMQMElOAr4MPCzJ6iQvG3VNU9gHeDHw9Pb38etJ/nDURU2yE/CF9u/6qzTnFD496J36jWZJUseegiSpYyhIkjqGgiSpYyhIkjqGgiSpYygsMkl+NvNSU677inaEzvVGjE3jPe28y5I8bgD7Ht+YI28mOSHJn7SPjxv1YGhJrp/LKLzt6/KnPdNLkrxnhnX+MslL2seHz2VU3SSfTPLgOax3cfs9i962w5O8d0O31a77iiQvncu66s9Q0Ib4H5ov0HxnUvuzgN3bf0cAHxhWQUnu9rfyq+rlVfWNjVHPCIwDXShU1aqqetV0K1TVB6vqI+3k4cAGhUKSPYHNqurbG7jeOM1YYrdvyHozOB6Y9vlqwxgKi1T76f7t7ZfeLk/ygrb9Hkne347f/ukkn1n3ibqqvlZV1/fZ3MHAR9pRHb8CbJNkpyn2e2w7PvxXkuyYZOsk17VDDpDkfu2n5nsmeXy77JeBI3u2cXiSTyQ5Azg7yVZJzktySftcDm6XW693kea+GG/qU9P5ae+dkOaeBZe0+z2vz7KbJXlHu5/Lkryybd8vydfa9uOTbNG2X5/kzT21Pbxt/6009774WpJ/ox3baLqakzwkybltbZck+R1gOfCU9stXr06yb/v/do9239v0bOva9jV/U7vdPwGWACvb9Q9McmrP8vsn6Tf8w2HAaT3L/SzJ29pewLlJ9mpf02+nGVRynWcBn2vXeWmSbyX5Is0XydZt64+SXNi+Lue29d4jyTVJxtpl7tE+l+3bb8dfn2SvPnVqDgyFxeu5NGO0P5rm0//b2zfy59J8+vxd4OU03+ydyWxH8Lwv8JV2fPgLgD9vhy0+HziwXeZQ4D+r6tfAh4FXVVW/Gp4ELK2qpwO/Ap5TVY8Dngb8c5J+A8hNq33T+Xfgj9san9dnsSNoxsR6bHtviJVJtqS5j8ALqup3acYU+6uedW5ta/sA8Hdt2zHAl6rqsTTDVTxoFiWuBN7X1vZk4GZgGfDfVfWYqnrnugWr6k6aN+7ntM/ticD1VXVLzzKfBFYBh7X34vgM8Ih1b77AS2n+DybbB7i4Z/q+wPlV9Xjgp8Bbae5B8hzgLT3LHQB8rv09e3O7nf1p7g+xzpeAvdvX5WTg9e1z+RhNGEHz+3ppVd3aTq8CnjLlq6YNYigsXr8HnFRVd7RvFF8EntC2f6Kq7mzHefrCLLY12xE8b6e5CRA0byrj7ePjaN6AaH9+OMn9gW2q6ott+0cnbeucqlo3bn+A/5tmOIBzaQJpx1nUPdnewAVVdR1Az/Z7PQP44Lp7Q7TLPAy4rqq+1S5zIvDUnnXWfdrufc5PpXmjo6rOBH44XWFJtgZ2rqpT23V+NYsxpD4OvKB9fGg7PaV2kMKPAi9qexhPov8Q5zsBEz3Tt9P2AIDLgS+2oX457fNNcx5hl/aQ0xNpQmSiPZTUW9cuwFlJLgdeB+zZth8PvKR9/GesH1Zr2MBDYJqaobB4TfVJeoM/YTP7ETx/3TM66h20o/RW1f8A40l+n+ZY9RVtHdONwfLznseHAWPA49tPvLcAW9IM0d37O77lDM9jpn1OtcxMr9lt7c/uObf67Wuqmufy//Jl4CHtJ/9D+E04TefDwIuAF9J8OFjbZ5lfsv5r2fv/eift820/4a97vk+h6QWsM9Xr/K/Ae9se11+s209V3QjckuTpNKHSG1ZbtjVpIzAUFq8LgBe0x8jHaD65XkTzh/vH7XHbHYF9Z7Gt04GXpLE38OOqunkD6/kIcBLtJ8Cq+hHw4yS/184/bKoVgfsDa6rq10meBvx2234LsEN7/H4LmmG7p/Nl4PeT7AaQZLs+y5wN/GXaE9ztMlfThNpD2mVeTNPzms4F655TkmcB205XczvW/+okh7TrbJHkPjSHa7but4P2jfpU4F9oRgP9fp/F1lu/qm6iCfR/oDkk1s9VwEOmmDeVA/jNG/mFwL7tc7wn6x+muz/w3fbx0knbOI6md3VKe5+BdR7KgIeTXkwMhcXrVOAy4FLg8zTHbr9HM778apo/sn+j+QP+MUCSV6UZQXQX4LIkx7Xb+gzwbeBammPyfz2HelbSvDGe1NP2UuB9aU40T/dJcCWwJMkqmjfaqwHaQxhvaZ/Dp5lh2OGqmqA5Z/CpJJfS/3DLcTT3676sXeZPq+pXba2faA973AnMNArsm4GnJrmEZmTbG2ZR84uBV7WHyf4fzY2CLgPWtiefX91nPx+n+eQ/1aGjE4APtiea7922rQRunOaKrDOZ3YeFXvvSBmX7geFNNCF8LnBJz3Jvonkd/xu4lfWdDmzFXc9z7NNuRxuBo6TqLpJsVVU/S/JbNL2HfQZ9H4n2SpiDq+rFg9yPZpbmOwNfq6oPTTH/3jTnmvaZ9Il9qu3tAvx7VT3rbta1BHhnVT2lp+2xwGv8vdl4DAXdRZo7PG0D3Av4p6o6YcD7+1eayxX/sOdkrUYgycU052v2r6rbplnumTSHpG4YUl3LaK7oOqyqvtTTvj9wzRSXSmsODAVJUuf1QMQAAAAiSURBVMdzCpKkjqEgSeoYCpKkjqEgSeoYCpKkjqEgSer8L4wqEPkw5cnxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot\n",
    "fig, ax = plt.subplots(1,1, figsize = (6,6))\n",
    "ax.hist(logK, bins = 40)\n",
    "ax.set_title('logK histogram')\n",
    "ax.set_ylabel('Count')\n",
    "ax.set_xlabel('log10 hydraulic conductivity (m/day)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Export the data for some kind of gridding\n",
    "\n",
    "df_Keep.to_csv(r\"C:\\Temp\\KeepRiver_logK.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "hydrogeol_utils",
   "language": "python",
   "name": "hydrogeol_utils"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
